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ABSTRACT 


This thesis is concerned with the one-electron 
approximation to the problem of finding solutions of the 
ochrodinger equation for the particles in a crystalline solid. 

The implications of the rotational and translational 
invariance of the Hamiltonian in the one-electron equation 
are reviewed. We see how the translational invariance leads 
to the Bloch form of the wave function and to the natural 
introduction of reciprocal space and Brillouin zones. We 
see also that the possible symmetries of the wave function 
are dictated by the rotational invariance of t!.e Hamiltonian. 

The anomalous skin effect, cyclotron resonance, and 
the de Haas van Alphen effect are reviewed briefly. We 
consider next some of the more standard methods for obtaining 
approximate solutions to the one-electron equation. These 
methods include the Wigner-Seitz cellular method, the Tight 
Binding procedure, and several variations of the Plane Wave 
methods. 

Three new procedures are then developed for obtaining 
approximate solutions to the one-electron equation. The 
first of these is a cellular method which provides conveniently 
E(k) for those special values of k which are highly symmetric 
in the first Brillouin zone. Next we develop a determinental 
procedure which produces an implicit approximation to the 
entire (E,k)relation. Finally, an interpolation procedure 
based upon the other two new methods is devised. 

The new methods are applied in some detail to the 
one dimensional empty lattice and the Kronig-Penny model, with 
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very satisfactory results in both cases. The methods were 
then applied to a plane hexagonal lattice with two distinct 
types of potential, and the results ere encouraging. Applica- 
tion of the new cellular method to a face centered cubic 
lattice proceeded until the capacity of the available computing 
facilities was exceeded. Some of the possible directions for 
future work are indicated. 
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CHAPTER I 


INTRODUCTION 


Since 1940 and earlier, physicists have been concerned 
with methods of obtaining approximations to the solutions of 
the Schrodinger equation for the particles in a crystalline 
solid. It is normally assumed that the atomic nuclei are at 
rest at fixed positions in the lattice, and there remains, 
then, an equation involving the wave functions of many electrons, 
the interaction of the ion cores upon these electrons, and the 
mutual interactions between the electrons. This many-electron 
equation has not been solved analytically, and a strictly 
numerical approach does not appear promising. 

The one-electron approximation to the problem assumes 
that each electron is adequately described by Schrodinger's 
equation with a potential energy term due in part to the fixed 
ion cores, and in part to some average potential due to the 
other electrons. The ion core portion of the potential that 
is used may or may not attempt to include the effect of the 
inner, tightly bound electrons. The one-electron approximation 
gives rise to two problems: the setting up of the appropriate 
Hamiltonian for the system, and then the solution of the 
resultant equation. 

Except for a brief discussion of the ton core portion 
of the potential in Chapter VII, we concern ourselves with 
methods of finding approximate solutions to Schrodinger's 
equation, given the potential. In Chapter II we see how the 
translational invariance of the Hamiltonian (assuming an 
infinite lattice) leads to the Bloch form of the wave function 


and to the natural introduction of reciprocal space and Brillouin 


zones. Following Mott and Jones (1958) we then see that the 
energy _E exhibits a band structure as a function of k where 
v= ek Ty. » and where ug (Pr) has the translational invariance 
of the potential. We then discuss the symmetry of the wave 
function which results from the rotational invariance of the 
Hamiltonian, and finally the symmetry of E(k) and the nature 
of surfaces of constant energy is investigated. 

To varying degrees of approximation, Chapter III is 
concerned with the theoretical bases for three phenomena which 
are more or less amenable to experiment, and which enable 
one to gain some knowledge of the geometrical characteristics 
of the Fermi surface of a metal. The three phenomena considered 
are the anomalous skin effect, cyclotron resonance, and the 
de Haas-van Alphen effect. 

In Chapter IV some of the more widely used methods 
for obtaining approximate solutions to the one-electron formula- 
tion of the problem are reviewed. The methods considered 
are the Wigner-Seitz cellular method and some of its variations, 
several Tight Binding procedures, and several variations of 
the Plane Wave methods. 

After considering a matrix equation which is equivalent 
to the Schrodinger equation for the problem, the standard 
variational procedure is discussed in Chapter V. Following 
this, three new methods for obtaining approximate solutions 
are developed: the method of Restricted Cellular Basis Functions, 
a determinental method, and an interpolation procedure based 
upon the previous two new methods. 

Chapter VI comprises several applications of the new 
methods which illustrate their use, and also their rates of 
convergence in some special cases. Chapter VII gives a brief 
review and comparison of some of the methods considered. 


CHAPTER IT 
THE THEORY OF ELECTRONIC ENERGY BANDS 


In this chapter we consider the nature of solutions 
to the Schrodinger equation 


Hy = Ey, Cea 
where H is of the form 

H= - y- + V(z) (2.2) 
and where V is periodic. That is 

V(r+An) = V(#) (2.3) 


where the i'th column vector of the matrix A is a; and the 
“elements of the column vector n are integers. For the sake 
of convenience we will treat r as a vector in three dimensional 
space, though all of the following material holds with but 
little modification in one and two dimensional spaces. Much 
of the material in this chapter is merely a reorganization 

of work which has appeared in many texts and papers. The 
main references have been Jones (1960), Mott and Jones (1958), 
Lomont (1959) and Wigner (1959). The notation, of course, 

has been made consistent. 

Since V° is independent of translations and orienta- 
tion, and since V(f) is periodic, the Hamiltonian H is 
invariant under two types of operation: translations by any 
integer linear combination of the basis vectors as of the 
lattice, and by any rotation and/or reflection which leaves 
the lattice (i.e. V) unchanged. 


Bloch Wave Functions 


Let us consider first the implications of the translational 


invariance of H. We denote by gq that translation operator 
which maps f(?) into f(f + An). That is 


gaf(r) = f(F+An). (2.4) 
Since 


g-e-f(r) = g-f(r+An) = f(r+A(ftm)) = g- f(r), 


(2.5) 


the operations (g=) clearly form a group G as m runs through 
all integer vectors. The group G is abelian, since 

Shim ~ Sm+A- touo) 
and H is invariant under the operations of this group; that 
is 


g,Hf(r) = He,-f(r). (2.7) 


Suppose now that 9(¥%) is a solution to (2.1) for the 
particular value of the eigenvalue E = Et. That is, we 


assume that 
Ho(#) = E'o(#). (2.8) 


If we now operate on both sides of this equation by the 
translation operator Zao we find 


g-Ho(P) = Hemo(r) = Ho(r+Am) ; (2.9) 
and 


goE'o(r) = E'o(F+Am). (2.10 


2a Le 
Ho(F+Am) = E'o(r+Am). Po) 


o(r+Am) then, is also a solution of (2.1), not necessarily 
distinct from g(r), corresponding to the eigenvalue E!, 
Then, given any solution 9(?) to (2.1), we can generate an 
infinity of solutions 9(r+Afi), for all integer vectors fi, 
all of which correspond to the same eigenvalue. Let us 
now order the vectors n in some fashion and label them Ny» 


and define the row vector 


oe) := |.--s0(e+an..),9(2raA),0(F+am), | 


(2,72) 


We now define the effect of operating on 9(r) with E, to 


g-6'(7) = [...,ec0(#+An_1),e,0(F+AR,) egp(PHAA,),...]. 


(2.13) 
But 


gro(F+An,) = o(F+A(K+45; )) 


o(r+An ;) (2.14) 
for some permissible vector nj. The vector gré is then a 


vector the components of which are merely a permutation of 
the components of 9 itself. We define 


gcd = Sz, (2.15) 


and the m'th component of $g is given by 
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[$,(#)] = g58,0(F). (2.16) 


Consider now a set of infinite matrices Gg, where 
each row of Gr contains exactly one 1, and each column 
contains exactly one 1, and the remainder of the elements 
of Gr are zero; that is, the (Gr) are a set of permutation 
matrices. Clearly, the elements of the column vector Geo(f) 
are again the elements of 9 in some permuted order; choose 
the Gr such that 


GeO = OF = Epo. (2,17) 


The set of matrices (Gp) obviously form a representation of 
the group G, and in fact, they form the so-called "Regular 
Representation" of G. 

Now, since (Gg) is a representation of the gyoup G, 


and since 
(sc-S"+)(sa.s-*) = SG=G5S"*, (2.18) 


then the set (sG-s~+) for any non-singular S also forms a 
representation of G. Because G is Abelian, it can be shown 
that S may be chosen such that every member of the representa- 
tion (86,874) of Gis diagonal. 

We had chosen Gg such that 


Ge6(F) “ 6-(F), (2.27) 
therefore 
(SGgS~")(S6(#)) = SG_(#). (2.19) 


Let us define the matrices 


=. 
J= = SGS (2.20) 


and the vectors 


Vz(F) = 86;(P). (2.21) 


We note that since the J's are diagonal, each set of numbers 
consisting of the n'th diagonal element of all of the matrices 
Jq» forms an irreducible representation of the group G. 

Now, (7) was chosen to be a solution to (2.1) 
corresponding to E=E'. If we say that ES is the unit element 
of G, then 


G5(F) = G(r), (2.22) 


and each component of 95 is also a solution corresponding 
to E*. and since the components of 9, are those of M5 in 
permuted order, each component of Or is a solution corres- 
ponding to E'. But the components of Vg are merely linear 
combinations of those of §-. In fact, from equation (o,21) 
we see that 


(Fe) in -) Snel@ide: (2.23) 


That is, the components of V_(P) are individually solutions 
to (2.1) corresponding to E = E!, 

Now, from equations (2.19), (2.20), and (2.21) we 
have 


Jnl s a Ve (2.24) 
and since each J is diagonal 


(Sp) mmlVsl., = Fel, (2.25) 


But operating on V5 with Jr is equivalent to operating on 
V5 with gg, since 
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= Vz. (2.26) 
Therefore, equation (2.25) implies that 


(Jz) ame Valin = SglFsl,,- (2.27) 


m 


Now, v=(#)]_, is a solution of Hy = Ey with E = E!; 
el¥s(#)),, = [¥s(F+ak)] 5 (2.28) 


and (Tp) sin is a number. Equation (2.27) then, implies the 
existence of a solution ¥(r) to HW = Ey with E = E!', which 
has the property that 


egv(r) = v(r+ak) =(Jz) wr). (2.29) 


(Je) am is the k'th element of one of the irreducible represen- 
tations of the group G of translations (g=). But the n'th 
element of the irreducible representation characterized by 

the triple of numbers 6 is given by 


e.= elB: jn, (2.30) 


where we now label the elements of G by a Single subscript. 
We note here that 


i (e1B+3n)(et8* dm) = elBs (IntIm) , 
(2 od) 


and also that for inequivalent representations we need 
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Equation (2.29) now yields 


ge, v(P) = ¥(FAR,) = etB-Iny(r). 


Consider now the function 


ug(#) = en 1B (AP) (5), 
ug(F+AJ) = en tB+(A™*(F4A9) Jy (443) 
= AB (AP) (5) 
= uz(r) 
If we now let 
(a-1)"8 = & 


and define 


ug(®) = ug(#), 


(2,90) 


(2.34) 


(2.35) 


(2.36) 


(2.37) 


equations (2.34) and (2.35) imply that we can write v(r) 


in the form 


w(z) = et Py, (3) 


(2.38) 


\O 


where ug (r+An) = ug (>) ; that is ug is periodic with the same 


period as V(r). This, of course, is the well known theorem 


due to Bloch (1928). (Jones (1960), Lomont (1959)). 
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clearly provide one acceptable volume V over which the vector 
B may range to generate all inequivalent irreducible represen- 
Pacions of the translation group G. It is aleo clear that V 
is not the only acceptable volume. In fact, any volume 

(or set of disconnected volumnes) which is equivalent to V 
will suffice, where the points with coordinates Bs and those 
with coordinates Bs + ern, are said to be equivalent. 


The Reciprocal Lattice; Brillouin Zones 


Let us set 
Nee ara = (2.39) 


and define the colum vectors of orBl to be the basis vectors 
of a lattice in "k-space" or "reciprocal space". From 
equation (2.36) we see that if 


B a ern, (2.40) 


Tieo 28, & JACLLCe DOING Of 4 BIMpLe cubic lartice of side 
ell, then 


(e)" = (f)7(278), (2.41) 
or 

k = (27B) 7A. (2.42) 
That is, k is a lattice point of the reciprocal lattice. As 


6 ranges through a unit cubic cell in B-space, k ranges through 
a corresponding rhombohedral unit cell in k-space. If in 
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B-space the unit cell has alattice point at its geometrical 
center, so does the corresponding unit cell in k-space. If 
the P-space unit cell has lattice points at its corners, so 
must the k-space unit cell. However, the most symmetric unit 
cell in B-space, a cube, does not correspond (except accident- 
ally) to the most symmetric unit cell in k-space. 

The most symmetric unit cell in k-space is easily 
found. The vector (2mB)!A, in k-space, defines a lattice 
point, and the equation 


-[(2rB)Ta J= (1/2) | (27B)*al® (2.43) 


wa 


defines a plane which bisects perpendicularly the line 
joining k = 0 to the lattice point specified by n. Simplify- 
ing equation (2.43) yields 


k-B'n = m-B!Bn; (2.44) 


the smallest volume enclosed by these planes, as n varies 
over all non-zero integer vectors, is the most symmetric 
unit cell in k-space. This cell is also, by definition, 
the first Brillouin Zone of the space lattice characterized 
by A. (Brillouin (1953), Jones (1960)). 

As k ranges over the first Brillouin zone in reciprocal 
space, B, in general, ranges over some unsymmetric unit cell 
in B-Space. Figures 1, 2 and 3 illustrate the relation 
among real space, reciprocal space and B-space lattices for 
the case of a plane hexagon. 

From equations (2.36) and (2.38) it is clear that E. 
may be regarded as a function of k or B. The advantage of 
considering E(f) lies in the simple volume over which it is 
necessary to allow 6 to range. It will be shown, however, 
that E(k) possesses certain symmetry properties, and these 
make it desirable to investigate the nature of E as a function 
of k rather than p. 
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FIGURE 1 


The Plane Hexagonal Lattice 


The hexagon shown is the Wigner-Seitz cell for the 


hexagonal lattice with basis vectors ay and an « 
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FIGURE 2 


Reciprocal Space Lattice 


The hexagon drawn in solid lines is the first Brillouin 
zone for the plane hexagonal lattice illustrated in 
Figure 1. f, and P, are the reciprocal lattice 
vectors. The point R and direction indicated corres- 
pond to the point R and direction on the polygon 

drawn in solid lines in Figure 3. The points P and Q 
(and directions) on the two rhombuses correspond to 

the points P and Q (and directions) on the two squares 


iy Bieure 3. 
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FIGURE 3 


The B-Space Lattice 


The three cells shown correspond to the three k-space 
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Energy Band Structure 


Assuming that V(r) is a bloch wave function, we have 
ik-r = 
vse *u(r) (2.45) 
where we suppress the subscript k in ug(r). Since u must be 


a periodic function of r with the same period as V, we may 


expand u in a Fourier series of the form 


ulB) + J) Cpe Prin Br (2.46) 
And we see that 
u(?+Am) = sy is Ck 
n 
= u(r), (2.47) 


since AB =I. Substituting (2.45) into Hy = Ey, we get 
V-u + 2ik-Vu + (E-k-k-V)u = 0. (2.48) 


If we now substitute (2.46) into (2.48) we get 


aaa (B-Ka+k--V) = 0, (2.49) 


where we have defined 


Zt 


= k-27B°n. (2.50) 


n 
If we assume an almost empty lattice, or 


ia oe Vly (2.51) 


we expect that 
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n(n )ee Cages Er (2.52) 
in which case 
eR ncGe ale p,emis Br cue Ks Po (2.53) 
and we expect, then, that 
E~ ([R|°. (2.54) 


The following material closely parallels that given by Mott 
and Jones (1958), the essential difference being that they 
consider only the case 8s = 0. In (2.49), then, we neglect 
all terms VC-, n # 8 as compared with VC=; we get 


7 ap -Prin'Br = be 
VCz = ab Cae (E-k- k=). (2.55) 


Multiplying both sides of (2.55) by e-72™°PY ana integrating 


over one unit cell of volume v, we find 


-, 2Tim-Br - = 
og ful#)e dt = C-(E-k- k= )v, (2.56) 


since the pe ansh ls being solutions of the Helmholtz equation 


with periodic boundary conditions, must form an orthogonal 
set for different values of m+B!Bm. The only terms which 
might survive, then are those for which n = + m. The product, 
when n = -m is ee Ti (em) Br which must be orthogonal to 
unity, and (2.56) follows. 

Let us now define 


v= = a/v) fu(z)e® Par, (257) 


Equation (2.56) becomes 
O-V. = c-(E-|-|*). (2.58) 
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If we now put m= 8 we find 
- 12 
ES |k5| + Ve (2.59) 


which agrees with our expected value E » [ee] 
Solving Tor 4 we find 


CaVq C-V- 


_ gees OYE a 
a E-|xz|° | &s| °-|R-| Sv; eee 


(2.60) 


This approximation for C= can only be valid if | | <4 | cs|, 
and this is clearly not the case when 


|*s|° - lal = 0. (2.61) 


Substituting equation (2.50) into (2.61) yields 
Pe TE ag = pe | ie 
k*B'(m-s) = 7 (m-s)-B B(m-s). (2.62) 


7% s = 6, equations (2762) is identical to equation (2.44), 
which defines the boundaries of the Brillouin zones. Taking 
s # 0 merely corresponds to translating the reciprocal 
lattice by an integer linear combination of its lattice 
vectors. 

We have found, then, that if J | ae NE the solution 
to (2.1) is closely approximated by 


y(B) = Coetks‘? (2.53) 
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and 


E= = les] “V5 (2.59) 
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except near the boundary of a Brillouin zone. 

If k is near the boundary of a zone, the above 
calculations are, of course, invalid. _-In this region we 
approximate to ¥ by the function 


W(?) = en (cue nr a8 age 
+C—e : (2.63) 


that is, we retain only the two large coefficients of u. 
Substituting (2.63) into Hy = Ey we find 


cuetKS-F(n_|z-|?-v)+c-01kii'?(z_|ic.|?-v) e Ge 


(2.64) 


If we now multiply this equation by e*a"r 


a unit cell we find 


and integrate over 


2 
c5(E+-| -V=)-C-V-_= = 0. (2.65) 
Similarly, for ean"? ae get 
2 2 
-CVa_= + 07 (E- |i -Vx) = 0. (2.66) 


To obtain non-trivial solutions C; and C= from (2.65) and 
(2.66), we find 
oe}: ei 2 
(E-|&-| °-v5)(E-|&s|°-v5) - \Va-s|~ = © 


(2.67) 


Near Ke = Ke, that is when k is near the boundary of a Brillouin 
zone, we find from (2.67) 
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Then, on the boundary of the Brillouin zone defined by 


ks = Ke, there is a gap in the permissible energies of width 


AE = 2|V- =| , (2.69) 


the nearest permissible energies being given by 
= 2 
E = | =| +V- +(v- -[ . (2.70) 


There is, however, no energy discontinuity across the zone 
boundary, since (2.68) is symmetric in m and 8. The energy 
discontinuity pointed out by Mott and Jones (1958) results 
from their considering s = 0 only. 

Figure 4 illustrates the situation in one demension, 
where V(x) = V(xtna), and V(-x) = V(x). We see that the 
(E,k) relationship for the empty lattice is given by the 
parabolas 
ent )* 


a ot = 0, $122) ews Cae s 


E = (k- 
A small periodic potential perturbs the (E,k) relationship 
in such a fashion that the continuous range of permissible 


energies 


2 
' ae ae <7 wt. 
T=2O0< E< 2 Geog y= SES za ee (2a) 


shrinks to the band structure 
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FIGURE 4 


Dispersion Relations in One Dimension 
The dashed curves are the parabolas E = (k-2n1/a)* 
and are the (E,k) relation for an empty lattice. 
The solid curves represent the perturbed (E,k) 
relation for a small periodic potential. 


20 


noltaxeqerd a 
vt gy 
= aie | 


: 918 BZavmMo Bedasb Cys Ye 
o L 
noitels (18) edd exe Dd Bm! 


it Jassstqe: gave Blog sal 
[lame &.9O% solve len ay 


“ al it se 
a ™~N 
= Y 
i ae! S 
ad 4 9|0 
~ i + 
ak, ZL ‘% 
~ 4 
™s Ps 
hy ra & 
Pres alo 
sain i. + 
- ie. 
— . 
a. \ 
a Ny # 
ee + 
a “i \ 
— 
™~ Pe 
~~ a 
—~ o“ O 
mM WN = 
Is = _~ 
- yy 
oa 
= NS 
* a a Flo 
~L \ | 
~sS e * 
— 
~ ~~ 
he 
pe u|o 
Ps sie, 
a = 
Aa NX 
aaa ie 4 & 
=— 7 
sell y~ oe hee 
a. Fa 
“7 
> a 
™~ a” 
~N 


rat | 


2 | 
¥<8 Se -1Ml > S+| Vy <2 iea hy Ue 
Wr? on® 
“3 cs Ivol< BE <. 2 = | v.| peas xe (2.73) 


We_note here, that in choosing wto bea Bloch wave 
function e“'Fuc(r), that it is not valid, even for very small 
periodic potentials, to regard k as a momentum vector except 
along those sheets of the (E,k) relation near E = Vale Only 
when s in equations (2.52), (2.53) and (2.54) is chosen to be 
O is this true. In general, 


BY = -iw = my-iet*’ wy, (2.74) 


and Vu# O unless s = 0. 
The free space Hamiltonian is, of course, a the 


empty lattice Hamiltonian is given by 
H = lim (-v°+ev). (2.75) 


The eigensolutions to the two problems are quite different; 


the (E,k) relationship for the empty lattice is E = |e | e Tor 


any integer vector s, whereas the free space solution admits 
only s = 0, Free Space is, of course, unable to exchange 
momentum with anelectron, am the momentum is indeed k. But 


a lattice, on the other hand no matter how small the potential, 
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possesses infinite inertia and is able, by means of Umklapp 
processes, to exchange momentum with an electron. The 
momentum of the entire system is, of course, conserved. 
(Peierls (1956). 


Wave Function Symmetry 

Since V(rtAn) = V(r) for any integer vector n, there 
exists a group of rotations, proper and improper, which leaves 
V and therefore H, invariant. 

We will superimpose a co-ordinate system on the lattice, 
and rotate the co-ordinate system leaving the lattice unchanged. 
Let r be any vector in the original coordinate system, and 
Let 


r = Sr (276) 


be the same vector in the rotated coordinate system. The 
matrix S, then, characterizes the rotation under consideration. 
Only certian matrices S leave V unchanged, and these might be 
determined in the following fashion. 

In the unrotated system the lattice basis vectors 
are as) which are the column vectors of the matrix A. Any 
point N given by the vector 


N - nya, = An (aa7%) 
. 
is a lattice point, where, of course, the components of the 
integer vector n are (n,). If we now rotate the coordinate 
system by the matrix S, the point N in the rotated system is 
given by 


Ni = SN = SAn. (2.78) 
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ee: | 
The basis vectors, called as in the new coordinate 


System are given’ by 


a, = Sa, Laa7o) 
and therefore 
At = SA (2.80) 
and 
Ee ae (2.81) 


Since S must be orthogonal. We find, then, that the point 
N is given by 


N =An (2.82) 


in the new coordinate system. If we now consider the vector 
N, and rotate it along with the co-ordinate system to get 

a new vector M, it defines a new point M the coordinates of 
which one found from 


Ta's = a'm (2.83) 
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in the rotated system. If we now insist that M also be a 
lattice point, for every integer choice of the vector n, then 
we are insisting that equation (2.83), rewritten as 


An = SA it (2.84) 
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have integer solutions m for every integer n. Equation 
(2.84) reduces to 


An = SAm (2.85) 


Any S which allows of such solutions characterizes a permissible 
rotation. 

In. practice of course, most if not all of the rotations 
which leave H invariant are readily determined by an inspection 
Se. gne Lattice. 


Let us consider, now, the point group of the lattice 


G, with elements EA? DS ee oe ears ak »N, where each gE, if a 
rotation and/or reflection which leaves H = VF +V invariant. 
That is 

g Hf = He f , n= 1,2,....... , N (2.86) 


By definition 
a f(r) _ 6(3 3) (2.87) 
Where SH is the matrix which characterizes the appropriate 
rotation am/ for reflection. As in thedevelopment of Bloch's 
1 


theorem, we assume that o(r) 489 a solution, for E=E , to 
Hy = EV; .1.e. 


Hp = E'g. (2.88) 
Operating by gE; we find 


g jHo(r) = He (>) = Ho(S ,r) (2.89) 
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and 
g Eto(r) z E'g .o(r) = Et9(S ,r) ; (2.90) 
therefore 
Ho(S .7) = Eto(S .r). [o.94) 
That is 


0 ,(2) = 6 ,0(F) 
is also a solution with E = E' to (2.1). Then, given a 
solution 9 we can generate a total of N solutions, not 
necessarily distinct, Ps» i = lye..N of Hv = By all corresponding 


te 2 = BY. 7 
Define, now, the row vector ® (r) to be 


Pie s - = 
6°(2) = [o,(#), 05 (F),...oy(#)], (2.92) 
; : —T 
and define the operation of Si. on Q' as 
6 -[go., go BOs (2.93) 
By? = LEY Pye By Paes + By Py 
We note that 


= 28.0 = £0 = 9 (2.94) 


o,(r) = g,o(r) (2.95) 
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so that ®,. is a vector the components of which are a permutation 
of those of ~. We now let Gs, y sole yaivan DN be a Bet of 
N x N permutation matrices forming the regular representation 


of G, and we chaose the subscripts so that 
@o(#) = 6,(%) = g, 0(z) (2.96) 


Again we perform a similarity transformation on the regular 
representation (G,) to get a completely reduced representation 
(J,), where 

J, = S8G,8™". (2.97) 


Again we note that G,° = §, (SG, 8~*) (86) = (S3,), that is 


IV =F, (2.98) 
where 
v= 3 (2.99) 
and 
V, = S9,. (2.100) 


It is clear that the components of V(r) are individually 
solutions to (2.1) with E = E!. 
Now, since the J's are completely reduced, that is, 


each ds is of the form 
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where the size of the block x) depends on k but not on m, 
and no further reduction is possible, the set ao) 
N forms an irreducible representation of G. This is true for 
each k, but the representations for different k may be 
equivalent. It can be shown in general that the regular 
representation reduces to give each t x t representation t 
times (Hamermesh (1962) ). Since G represents the point group 
of a lattice, no irreducible representations are of order more 
than three. 

Let us assume, now, that g, is the identity operator. 
One of the irreducible representations, say the first, must 


be the identity representation; that is 


g@) er ee ree (2.101) 
And since 

3.4, (#) = ¢9,(#) = 9,(5,2), (2.102) 
then 

s]) [9,(2)1, = (W(s,7) 1,. (2.203) 
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But 3) = 1 for all m; equation (B°9077., then implies that 
the first component of V,(r), call it @(r), is a solution to 
(2.1) with E = E' such that 


a(S_r) = O(r). (2.104) 
Sw? here is any rotation or reflection which leaves H invariant. 


Another irreducible representation, say the second, 
might be an alternating representation, that is 


g(2) = 1 for some m, l<mM 
g(2) = -l for the remaining m. (2.105) 


Then the second component of v,(), call it e(r) is a solution 
such that 


(2.106) 
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If the third irreducible representation is also l by 
1, the 33) must be of the form 


363) 2 et %m (2.107) 


for some set of real numbers (a). In this case lv, (¥], is 
a solution which has the property that 


[v1 (s,7)], = ev m ly, (F)J, (2.108) 
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Let us assume, now, that the (s)) form a 2x2 
irreducible representation Rename the fourth and fifth 
components of V,(r) é(r) and ¢(r) respectively. Then 


ay fe) fs yak + oy r08 
a = 1) ; ; (2.109) 
ef Ma ent + (t)) 05, 


but this must be the same as 


A(z) c(S_(r))}. (2.110) 


That is € and C must be pair of solutions such that 


e(s 2) = (TS) yay e(3) + (a4) ye (8) 


c(s@) = (Tt) )or €(3) + (900 e(R). (2.222) 


Thus, the irreducible representations of the point 
group of the lattice determine the possible symmetries of the 
solutions to Hy = Hy. 

The above considerations refer to solutions of 
Hy = Ey which are not necessarily in the Bloch form. If we 
do insist that our solutions be of the form Sate then 
we are interested only in those rotations which transform 
solutions into themselves, or linear combinations of other 
solutions belonging to the same k. That is, we must restrict 
the group of operations to those which leave k invariant, or 


change it into an equivalent k 
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If we rotate co-ordinates by a rotation characterized 
by 5S, the new coordinates of any point defined by r are given 
by 


r = Sr (2.76) 
And since 
Awa’ = SA, (2.80) 
then 
BB = Bs! . (2.119) 
pance.1t-is 
(a. \" = ea" (2.113) 


which characterizes the reciprocal lattice in the new 
coordinate system, it is clear that 


K,, #2 Sk. (2.114) 


We insist, then, that given k, the rotations Sn be confined 
to a subgroup G(k) of G, such that 


zi 


Sk = ktonBon. (2.115) 
Clearly, for most vectors k in the first Brillouin zone, the 
subgroup G(k) reduces to the identity element. However, for 


vectors k along lines or planes of symmetry inside the zone, 
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there are other elements of G which leave k invariant. And 
if k is also on the zone boundary, there will be other elements 
of G which transform k into an equivalent vector in reciprocal 
space. 

Again, it is not too difficult to determine the sub- 
group G(k) given k, by an inspection of the reciprocal lattice. 


The S etry of E(k); Surfaces of Constant Ener 
If wy is of the Bloch form, then Hy = Ey yields 


V-uteik-vut(E(k)-k-k-V)u = 0. (2.48) 


If we now perform on the space coordinates a rotation 
characterized by the matrix S, we are in effect making the 
change of variable 


af (2.76) 
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As we have seen, the space rotation induces a rotation in 
reciprocal space also characterized by S, so that 


t=c.e 4 (2.114) 


Since equation (2.76) implies that 


oe (2.116) 


the various terms in (2.48) undergo the following transformations; 


v . viv = (y')*ss*y' Za ( 147 i (y')e 
(2,117) 
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the well known results that neither vy nor the length of a 
vecotr are affected by rotations, proper or improper. Nor is 
the term k-V: 


(2,119) 


If we now insist that S be one of the rotations which leaves 
H invariant, equation (2.48), after dropping the primes, 


becomes 


V°u(sS'r)+21k-vu(s'2)+(E(S kz) -k-k-V(F))u(Stz) = 0. 


(2.730) 


This equation is of exactly the same form as equation (2.48), 
and since periodic boundary conditions are unaffected by 
rotations, we see that symmetrized solutions u(r) to equation 
(2.48) transform according to the irreducible representations 
of the entire point group of the lattice, and that E(k) is 
invariant under all operations of this group. 

The symmetry of E with respect to i and the fact 
that E is an analytic function of k except where different 
energy bands touch in k-space (Jones, 1960), greatly simplify 
the determination of the (E,k) relationship within the first 
Brillouin zone. In the case of the plane hexagon, for example 
one need consider only those k's within the triangle bounded 
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The energy over the remainder of the zone is easily generated 
by suitable rotations and reflections, given the energy for 

k inside this triangle. In general, a knowledge of E(k) where 
k varies over some small fraction of the Brillouin zone 
provides one with E(k) over the entire zone. 

Further, since E(k) is analytic and invariant under 
the operations of the point group of the lattice, surfaces of 
constant energy must intersect at right angles the various 
planes of symmetry of the Brillouin zone. And on the boundaries 
of the zone, non-degenerate energy bands must intersect the 
bounding plane at right angles. 

The solid curves of Figure 4 provide an illustration 
in one dimension. The only "plane" of symmetry within the 
first Brillouin zone is the point k = O, and the boundaries 
of the zone are at k = t1/a. 
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CHAPTER III 


EXPERIMENTAL METHODS FOR DETERMINING FERMI SURFACES 


In this chapter we consider three phenomena which 
are more or less amenable to experiment, and which enable one 
to gain a knowledge of some of the geometrical characteristics 
of the Fermi surface of the metal. To varying degrees of 
approximation, we concern ourselves below with the theoretrical 
bases for the anomalous skin effect, cyclotron resonance, 


and the de Haas-van Alphen effect in that order. 


The Anomalous Skin Effect 

Consider the conduction of electricity near the 
surface of a semi-infinite slab of an isotropic metal. The 
surface of the metal is the x-y plane, and the current is 
being induced by an electric field of frequency w, plane 
polarixed.in the x direction, normal to the surface of the 
metal. Assuming that electrons which collide with the surface 
of the metal are reflected with kK, and Ky unchanged, and kK 


replaced by “ko: Pippard (1960) has shown that the surface 


impedance Z in the x direction is given by 


co 


Z = 8iw (Sea) 
ey +4riwo(v) 


Here, the conductivity is defined by 


J =oal(v)é (3.2) 
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and EY and Fr are the Fourier components of € (z) and the 
current density J(z) (both in the x direction); that is 


G(z) = [ & e*¥av (343) 


and 
J(z) -fa ar Eis (3.4) 


If the free path # of an electron is short enough, then the 
only values of v which make a significant contribution to Z 
in (3.1) are those for which 


OF lee ae (4.5) 


In this case o(v) may be taken to be constant over the range 
with a value equal to the D. C. conductivity o, and equation 
‘oe is easily integrated to yield the normal skin effect 


Z = (1+i)N(27w/c) (3.6) 


n 


If, however, the free path of an electron is long enough, the 
significant values of v in (3.1) satisfy 


ny an (3.7) 


Pippard (1960) has shown that as v , o(v) is given by 


o(v) = a prose, (3.8) 
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where ROK) and the range of integration are determined as 
follows. Consider the curve formed by the intersection of the 
Fermi surface E(k) = E,, and a plane normal to k.. RUKL) ig 
the sum of the absolute values of the radii of curvature of 
all points on this curve the tangents of which are parallel 
to the kK axis. The minimum and maximum values of ky are; 

of course, given by the extrmities of the Fermi surface in 
the ky direction. In this case where vi—7™ the extreme 


anomalous limit), Z is given by 


oO 
‘tal 
Lug = Bie | av me. (3.9) 
Oo ¥ +ia~ 
where 
ror] — se R(k A (2 4 . * 
_— ab } K e / @ ae / 
TH [ / y 


Equation (3.9) is easily integrated to give 


i Sirw cf oy 2.937) 
Pixs ~ a3) Re (1+iN3) e (3.11) 
oD © 


It was shown by Reuter 2nd Sondheimer (1948) that 
assuming completely diffuse scattering at the metal surface, 
as opposed to specular reflection as above, leads to 
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Z = in sul [in(1+4nien(v) av], (3.12) 


which, in the extreme anomalous limit yields 
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The impedences given by equations (3.11) and (3.13) differ by 
a factor or 9/8 ; the second seems more plausible from a 
consideration of the metal surface. In either case, Z depends 
only on the frequency of the applied electric field, and 
geometry of the Fermi surface. Measurement of Z for a variety 
of crystals orientations provides one with,, and hence 
fixe: )dk 
J o 

As implied by Chambers (1956), reconstruction of the 

Fermi surface from a knowledge of this integral would indeed 


k 


be a formidable task, and is not even possible unless the 
Fermi surface is both convex and centro-symmetric. The 
anomalous skin effect, then, would appear to be a more suitable 
tool for checking model Fermi surfaces derived by other means, 


than actually constructing the Fermi surface of a metal. 


Cyclotron Resonance 
Let us first consider the phenomenon for the case of 


a classical electron. If the electron is subjected to an 
electromagnetic field of electric component & and magnetic 
components 4, then the force on the electron is given by 

1 xv) (3.14) 
Cc 

where v is the electron velocity. The total energy E, then, 
absorbed by the electron is given by 


E(t) = P(e) -H()at (3.15) 
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since vAxv is zero. 
The equations of motion of the electron are given by 


~ % 


setting the right hand side of (3.14) equal to mv. If &is 
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circularly polarized in the x-y plane, th 
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Ey 
f= 0 (3.17) 


and if Wis of magnitude ¥ in the z direction only, then the 
eugations of motion of the electron are 


t+ Ov. = g  coslat+9) 
7 oO 


v- QV. = _o sinfuwt+e) | 


Here, the "cyclotron frequency"? is given by 


Qo = e M/me (3el dy 
and 
o = eb fut?) /m. (3.2 
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HBr Es O, equations (3.18) are homogeneous, and their 
solutions are 


= AcosQht + B singt 
Vy = -B cosnt + A sinnt 
Vv, = constant (3.21) 


where A and B are arbitrary, and the electron describes a 
circular orbit with frequency 2 in v-space. If E+ 0, 
particular solutions of (3.18) are given by 


v. =o gin{wtt+? )- o 7.4, \ 
Xx a TS ( ; a0 [ sin/at+9) + gsin(Qt-9) ] 
w-2 w= wt+Q 
v. = _-¢_ cos(wt+9 ) + g [cos(att+9)+ cos(at-o ) ]. 
Seon 29 wH0 wD 


(3.22) 


Tf wf O, equations (3.15), (3.21) and (3.22) yield 


E(t) = eo [i-cos(w-a)t + cos20 -cos!(w-a)t+29) | 
20 Cm-Q)e MceQe 


+ Aesin((w-Q)tt+o) - Belcos@ -cos((w-2)tt+g) J, 
w= w-2 
(3.23) 


and we see that the total energy absorbed by the electron is 
a periodic function of time. 
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If, on the other hand, w=, equation (3.23) is 
replaced by 


2 
ear a 
B(t) = bq Lp + Sag t] + e(acose-Bsine)t, (3.24) 
: 3.24) 


and in this case the energy absorbed by the electron is a 
quadratically increasing function of time. 

Let us now assume that for some effective mass m* 
(Chambers'(1956) "orbit mass"; Pippard's (1960)"cyclotror 
mass") that we can rewrite equation (3.14) in the form 


Ak = e(xv)/c (3.28) 
where £ is now zero and v is given by 
7 = V-ER(K (3.26) 
Ay -E(k). 


In rectangular coordinates (3.25) and (3.26) give 


mk, = - oh, 82 
- he 9k. 
y 
ae eH , OF 
th, = * fo | Bk (3427) 


and k, =O. If curvilinear coordinates k, and k, are set up 

in k space such that kK, 14s tangential and kK, normal to the 
trajectory described by equations (3.27) then these equations 
imply that 


Ak a rr | (3.28) 
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since on /k, = ©, where we have taken oe ak, 3/8 rk) to 
have magnitude unity. 
. a \ 1 » , 
Equation (3.28) may be rewritten as 
fe qk 3 
va  . AL os 4 ic 
' = a ~ a. 1 \ 4+ t 
d Gs ee § d= vw ( 3 od / 


and therefore if the trajectory in k-space 


orbit time T is given by 


nee J Kn dik, . a. 
= ee wi O8ar 
But 
f (dk, dae, = dA Pe 


where A is the area enclosed {in k-space) by the orbit, and 
SO 


r= oc dA. (3.32) 
ex dE 


The cyclotron frequency Q, may then be written as 


Q= 27 = eX (3.33) 
<7 aa 
Ki "4 m°’c 
where 
A~ A e -~ 8 
nt = 2 oA (3.34) 
Tr as 


If, now, we assume that for smell & we can neglect the 
quantization of E and A and write 


oe 
wai 
at Ca. we HVS: devet ave. ew. “rad a = is. z 
ies. 7 Minn 
om ie Viaoit ook 2 aaah 


c 


y 


é b= eal OF he ‘3 ag er. 4 Py *Troynst =a Siit- tr 


- ys covlg-et F 


ba .Jhdge sfi7 yd Se ee mb) eeCLS 


Ee; | % AE a) 
oH ke 
- oe 


ce 1S724N sd) fend een Gi 


unt #&) 


8g tuetmen? aso s 


Ho 


E = A= /om (3435) 
and 
A = 1k '3..36) 
for an otherwise free electron, then 
di = 7m (S237 
eka) Ane ~ 
and 
m= mi 3.38) 


as expected. 
If, however, E and A are more complicated functions 


of k, then measurements of 2, the frequency at which resonant 


a 


absorbtion of energy occurs upon application of a circularly 
polarized electric field, provides m* from equation [3.32), 


and thus dA/dE with the aid of (3.34). 
Now, m* depends on both E and ko but since only tho-- 
electrons near the Fermi surface are able to absorb energy, 
Aes 
described, however, provides no means of singling out a 
particular k,, and so the measured m* will be an average value 


any m* measured will be approximately m* (Esk ). The procedure 


over k, of m* (Ep, ue 

Chambers (1956) has pointed out that if the applied 
magnetic field is almost parallel to the surface of the metal, 
then only those electrons for which * in the direction of 
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the field is zero, or very small, will return to the 
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frame 


place in the surface layer and thus be able to resonate. 


this fashion one should be able to measure m* for a 


localized 


~ 


group of electrons, rather than obtaining an average of m* 


for the whole Fermi surface. 


The de Haas-van Alphen Effect 


During experiements concerning the magnetic susceptibi- 


lity of single bismuth crystals at low temperatures, 


de 


and van Alphen (1930) discovered that the susceptibility 


Hooe 
aas 


varied in an unusual oscillatory fashion with the applied 


magnetic field. Following the treatments given by 
(1957) and Chambers (1956) the origin of the effect 


determined in the manner outlined below. 


Given an electron traveling on a curve on a surface 


Shoenberg 


may be 


of constant energy in k-space, the applicatim of a constant 


magnetic field will alter the orzit in k-space to allow only 


those orbits corresponding to the permissible energy 
E,- It has been shown by Onsager (1952) that if the 
in the direction, then the allowed orbits in k-space 
for which A, the area of the projection of the orbit 
Ky - Ky plane, is given by 


A = (nt+y,2rek/Ac 
where n is an integer and yY is an undetermined phase 
The value Y = 1/2 gives the correct energy for an em 


For an infinite empty lattice, then, the Fer 


or 
upon application of a magnetic field in the direct 


with k. as their axes. Each of the quantized states 
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of sufficient degeneracy that the total number of states 
previously available in the sphere is now available on the 
cylinders (Pippard {1960)}. For a non-empty lattice the 
Fermi volume will be broken up not into circular cylinders, 
but into tubes of variable cross-section which may be inclined 
to the k, direction. 


The total energy W, of the electrons in such 


situation is given by 
oO 
. ereK) py, ) [2vdkz J 
W - fy (2s ] ey Ae Ky [ Pe 63 
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Here v is the volume of the metal, E. the energy of the n'th 


n 
permitted state given and k,, and ereH#/nce, the area between 
3 


permitted orbits, times avdk,/(2r) gives the number of states 


Sth, Ke The upper limit of the summation need only be taken 


. b ir f 4 \ ay a 2 4 r 
as the largest n for which EB NALk, “Ep, the Fermi energy. n, 


ear) 


of course, depends on Ko and we assume here that FE, i 


ia] 


independent of ®, which is approximately true except for very 
high values of KH. 
Thus, as 4 increases, a term drops out of the inte- 
\ 
} e 


grand of equation (3.40) whenever E (H,k,) grows to E., or in 
tz oo TT ee 
other words, whenever A \E yo HK) grows to Aves I,K... Let us 


assume now thatK = ! makes Aj equal to A(E5,0,k, ). From 
equation (3.39) 
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Let us also assume that increasing *® from | to bt! makes 
Ay 4 equal to A(E p30, k ops Then 


since the left hand sides of (3.41) and (3.42) are both equal 
to A{E,,0,k,), we find that 
= = ere (4,43) 
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The integrand of equation (3.40), then has discontinuities 


spaced 4 apart in wae that is, the integrand is an oscillatory 


function of 1/4, with "period" 2re/AcA{E.,0,k,). 


if Z 


Now, as } increases, the “tubes” in k-space containing 
the permissible orbits will expand eccording to equation (2.25). 


Eh 


it. is clear that the expansion most of the tubes can vary 
W only smoothly. There are only two conditi 
expansion of a tube can cause a discontinuity in the integrand 
of W: those tubes for which the permissible range ky 
(determined by Ee ye so small that any increase in (and 
therefore A) causes it to vanish, and those tubes which are 
tangent to the Fermi surface at some point, so that any increase 
in A causes the tube to split into an upper and lower tube. 
For convex Fermi surfaces only tne first condition applies and 
therefore only that value of k, for which A{E., ae is a 
maximum will cause discomtinuities in the integrand of W as 
ich aremt convex, or even 


mr 


kW increases. For Fermi surfaces wh 
of one sheet, all the values of k, for which A(E,,0,k, ) is 
eer err in the 


09 


maximum or a minimum will cause 
integrand of equation (3.40). (See figure 5). 
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FIGURE 5 


The de Haas-van Alphen Effect 


The two cylinders at the local maxima of the cross 
sectional area will disappear with increasing 

magnetic field. The minimum in the cross sectional 
area of the Fermi surface has just caused one cylinder 
to break into two. 
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If we now c 


w 


ll one of these extremal cross-sectionel 
t 
areas of the Fermi surface A , we see that W, and therefore 


Oo 2 en en 4 ae) » : ‘ 
xX = - 8°W/8X", will oscillate with period A in 1/4 given by 


= sof t f S, 
O& = oreyhcaA . 3.44) 
One must point out that if several cross-sectional extremes 


of the Fermi surface exist for a given crystal orientation, 
then x will oscillate with several periods, and a complicated 


beat structure in y(1/4) may result. 


We note here that the difference AH between os and 
z 
His 33 when A is near A , is approximately given by 
A Os aed AoA 
Ak = GE(A ) ak = ore GELA = e Hh , 
ae < , t 
A Ae dA m (A Je 


where m* is the “orbit mass” encountered in cyclotron resonance, 
given by equation (3.34). And since 2, the cyclotron frequency 


is given by 


then 


E 


that is, the cyclotron frequency 2 is simply that frequency 
which provides the appropriate amount of energy A& for 
transitions between adjacent quantized orbits. 

In the above analysis we have implicitly assumed that 
‘the temperature T was zero. if T is not zero, but large 
enough such that kT~A4E, the energy difference between successive 
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quantized orbits, then, as the cross-sectional area A of an 
orbit increases with through a somewhat fuzzy Fenmi 
surface, the change in the integrand of W (from equation 
(3.40)) will no longer be discontinuous, but will still be 
marked if T is small enough. The period, then, of W will 

be unaffected, but the amplitude of the oscillations will be 
reduced. One :expects, then, that the amplitude of the 
oscillations of y (which is equal to 8° R/ay4o Pee es 6G, 
where F is the free energy) to depend on 


! 
kT = kTm*(A_)c (3.47) 
AE en 
when “TA O. Lifshtiz and Kosevich (1954,1955) have shown 
that this indeed is the case; they find for tlhe oscillatory 


term in F 


cists | T#2/2cos[Aca' /eh) +6] exp-[ 2r?kTm* (A')c/enn] 
(3.48) 


Thus for various orientations of the magnetic field one can 
find from measurements of the period in 1/Hof xy, the extreme 
cross-sections A’ of the Fermi surface, and from amplitude 
measurements one can find m*(A'), that is da‘ /dE. We note 
that m*(A') from equation (3.48) is the effective mass for a 
particular k-space orbit, rather than an average value from 


several orbits. 
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CHAPTER IV 
STANDARD METHODS FOR CALCULATING ELECTRONIC ENERGY BANDS 


In this chapter we examine briefly some of the more 
widely used procedures for determining (E,k) relationships. 
The methods considered can be classified roughly into three 
groups: cellular methods, tight binding methods, and plane 
wave methods; and this is the order in which they have been 
presented. The classification, however, is far from unique, 
Since many of the procedures contain aspects from more than 
one of the main groups. 

In addition to individual references to them, exten- 
sive use has been made of the review articles written by Reitz 
(1955) and Pincherle (1960). 


The Wigner-Seitz Cellular Method 
The cellular procedure described below was first 


suggested by Wigner and Seitz (1933, 1934). With this method 
the crystal is divided up into the most symmetric unit cells, 
commonly referred to as Wigner-Seitz cells, and inside each 
of these cells the crystal potential is approximated to by 

a potential which is spherically symmetric. 

Using this approximate potential the Schrodinger 
equation is now separable, and the solutions to the angular 
parts are, of course, the spherical harmonics wy Pst) The 
radial solution f,(r,Eg) satisfies the differential equation 


r + (Ep - V(r) - 4(44+1)/r°)f, = 0, (el) 
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and so 


¥-(2) = » agi, (9.9)2, (r,B-)/r. (4.2) 


Lym 


The solution (4.2) is, of course, valid only in the cell 
centered at r = 0, but similar solutions ¥7(r = R, ) are valid 
in cells centered about the lattice points given by R,. 

In principle, the procedure would be to determine 
the boundary conditions that ¥ must satisfy on the surface 
of a Wigner-Seitz cell by means of the Bloch condition 


Wc(FHR,) = eA My (3); (4.3) 


k 
these boundary conditions wuld then be used to determine the 


¢m in (4.2), and, as usual, one would find that only for a 

discrete set of E-'s, given k, could one find non-trivial 

solutions for the Bm 
The boundary conditions are easily determined as 

follows. Let R be a lattice vector corresponding to the 

center of a nearest neighbour cell to the cell under consideration 

(centered at r = 0). Let ry define a point on the face between 

these two cells, and let LP, define a corresponding point on 

the opposite face of the zero centered cell such that 

Py = ry - RR. Let € be an infinitesimal vector in the direction 


of FR so that r is inside the zero centered cell, while 


ac 
ry + € is inside the neihbouring cell. For continuity of ¥ 
we require that 


WclB,-<) = We(F, +e), (4.4) 
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but rte = r,.+R+e, and so 


As 2 
Vg(Py-2) = Vg(Fotfte) = eM RY (B42), (4.5) 
and as 0 we find 
e(#,) = Ry (z,) (4.6) 


Similarly, for continuity of w we require that 


ve(Fy-e) avg (Fyt8) ask av (Fate) 
an, an, 7 an, 


hwere ny and No respectively denote the outward normals on 
the mutual face of the zero centered and neighbouring cells. 
As esO we obtain 
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In principle the boundary conditions (4.6) and (4.8) 
Should be applied to ¥ to determine the E- and the coefficients 
a, in the expansion (iia). In practice ois computational 
difficulties arising are prohibitive, and Wigner and Seitz 
(1933, 1934) made the following approximations. 

They used this procedure for k = 0, and retained 
only the first term in the series expansion (4.2). The wave 


function, then is given by 


Vel) = ¥e(r) = Flr) (4.9) 
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where F(r) = a f (r,E-)/r. 
©°O oOo kK 

The boundary condition (4.6) is now automatically 
satisfied, since k = 0 and v is a function only of the length 
of the radius vector, r. The Wigner-Seitz cell was then 
replaced by a sphere of radius ry and volume v, equal in 
volume to the unit polyhedron, and the remaining boundary 
condition applied to this sphere. In this case equation (4.8) 


implies that 
aF(r,) _ o, (4.10) 


For non-zero values of k, Bardeen (1938) constructed 
a procedure to find the eigenvalues of other states by means 
of perturbation theory. The wave function is written in the 
Bloch form Paarl ca and substituted into Schrodinger's 
equation to get the differential equation satisfied by u(r). 
The term in kYu in this equation is then treated as a perturbation. 
The method involves expanding u(r) and E as power series in k, 


and it is shown that 


E- = Et alz|® + o( fx"), (haa) 

where 
a = y[(r/P)(aP/ar)-1)],, , (4.12) 
Y= (47/3) rel v(r)]*, (4.13) 


and where P is a solution vanishing at the origin to 


p''+ (E,-V-2/r°)P = 0. (4.14) 
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Equation(4.11) is of course valid only when the 
expression kVu is small, and this is the case, in general, 
only for small | ic| 


The Expansion of y in Lattice Harmonics 

This method, developed by Slater (1934a, 1934b) and 
Von der Lage and Bethe (1947) is essentially a modification 
of the Wigner-Seitz cellular method. 

The crystal potential inside one Wigner-Seitz cell 
is assumed to be spherically symmetric, in which case equations 
1) end (432) apply. 


fy + (Bg -V(r) - e(st)/r®)e,=0, (4.2) 
vel?) = } am VYO.9)£,(2,Ez)/r. (4.2) 
4,m 


It was suggested by Slater that instead of approximating 
to the Wigner-Seitz cell by a sphere of equal volume in order 
Lo permit the convenient application of the boundary conditions, 
that instead one limit the number of 4's used in the expansiorm 
(4.2), and then apply the boundary conditions at a number of 
specific points on the surface of the Wigner-Seitz cell itself. 
It was found, however, that the resulting values E(k) 
were quite sensitive to the positions on the unit cell at 
which the boundary conditions were applied, unless the maximum 
fused was quite large. For this reason only points k posses- 
sing a high degree of symmetry in the first Brillouin zone 
should be considered with this method, since for these values 
of k many of the spherical harmonics are of unsuitable symmetry. 
This, of course, implies that their coefficients must be zero, 
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thus permitting one to use a fairly large 4 and still have 
only a reasonable number of terms in the expansion (4.2). 

A further improvement was made by Von der Lage and 
Bethe (1947), by rewriting the expansion (4.2) in the form 


ae =) Pyke 4(O00)2, (eBg)/r. (4.15) 


4,m 


Here K. 446859) is a "lattice harmonic", or Kubic harmonic" 
for lattices having cubic symmetry. Kp 4 (8.9) is a linear 
combination of the spherical harmonics (9 ,®) which trans- 
forms according to the irreducible representation labeled 
s of the appropriate subgroup (dictated by k) of the point 
group of the lattice. There may of course be more than one 
such linear combination for a given #, hence the subscript t. 
The lattice harmonics are independent of the crystal potential 
and to some extent, of the lattice itself, and therefore 
having once been found, they can be used with a variety of 
potentials in crystals having the same symmetry properties. 

Von der Lage and Bethe tabulated Kubic harmonics for 
the full cubic group. These apply for k = 7/a(1,0,0) in body 
centered cubic lattices, for k = 7/a(1,1,1) in simple cubic 
lattices, and k = O in cubic lattices, simple, body centered, 
or face centered. Howarth and Jones (1952) have tabulated 
kubiec harmonics for other symmetry points, and Bell (1954) 
has given lattice harmonics for various structures. 

The use of lattice harmonics reduces considerably 
the number of terms for a given 4 in the expansion of Y. 
The boundary conditions are now applied at selected points on 
opposite faces of the Wigner-Seitz cell to determine which 
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values of Ee permit non-trivial solutions for the coefficients 


Die: 


Altman (1958) has given 


applying the boundary conditions. 


enough points on the boundary of 
determine the b 
N, 


te? the procedure 


of points on the boundary of the cell, 


an alternate procedure for 
instead of picking just 

the Wigner-Seitz cell to 

is to pick some large number 


and then minimize 


the sum of the squares of the errors for each boundary condition, 
tat 285 The 


procedure then reduces to the solution of a matrix eigenvalue 


for each point, a least squares technique. 


problem. 


Linear Combination of Atomic Orbitals; Tight Binding 
The following method was first proposed by Bloch 


(1928) and comprises an attempt to approximate to a crystal 
wave function by a linear combination of atomic orbitals. 
It is assumed that the atomic orbital ee) satisfies 
the equation 
- E° 4.16) 
-V Xn + WX, a EX? (4.16) 
and also that we can write the crystal basis function AT) 


in the form 


Ven(®) =) eFax, (2) (4.17) 
J 
where 
X46) = Xy(P-Hy? (4.18) 
The R. are, of course, lattice vectors. 
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We note here that Vin is in the Bloch form, since we 


may rewrite equation (4.17) in the form 
Ven(B) = eM PY eA MRS-F)y (3) (44.19) 


where the factor involving the sum over j clearly has the 


period of the lattice for an infinite crystal. 


We now assume that the crystal wave function V(r) 


1 given by 


¥e(B) =) Avg, CF) (4.20) 


n 
where the sum over n is to be taken over the number of 


degenerate states of the atomic level Eo Substituting (44,20) 


into 


VW; + V(r)¥e 7 EVs (4.21) 
where V(r) is the crystal potential, we get 
z _ 
) Ag(y Vo, + (E-V)¥g,) = 0- (4.22) 
n 
But 
fice 
V Ven a i. IV Xn 4) (4.23) 
which, using equation (4.16) becomes 
2 _ ik-R _We 
vg, = Yate (wma | (4.24) 
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where 


W(p-R,). (4.25) 


W.(r 
Ae 
Then, substituting (4.17) and (4.24) into equation (4.22) we 
get 


ik-R. ; o - 
> Aral de 3 (W5-V+Ep -E) Xn 3 ) = 0. 
n a) 


(4.26) 

We will now multiply this result by Ras and integrate over 
all space. We assume here that the Xn 4? Pitts 2 oS gee ave we form an 
orthonormal set, and also that the "overlap" integrals may be 


neglected; this is, for this calculation we assume that 


x ba lo 
Tee a ade (4.27) 


Equation (4.26) now yields 


Destin [(Ee - E°)é YB RR) [Ca Sh 87] 


n’~mn 
n J 


Since we are summing over j, the left side of (4.28) is 
independent of i; we will drop the subscript, then, to get 


» Ay [(E, “E)6 sis os LER ifx (W -V) )Xn jor] = : 
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4.29) 


For a non-degenerate s level the sum over n contains only 
one term; only the value m = n need be used then, and we find 
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that 


° ik-R, * 
Ee = Be - ) et j ice (W5-V)X, 597 (4.30) 


as an explicit dispersion relationship. 

On the other hand, for a triply degenerate p state 
the sum over n contains three terms n = ny igy Ag: we must 
therefore also use equation (4.29) with m = Ny »M5 Nz. The 
condition in this case for non-trivial solutions A, is the 
vanishing of the determinant of coefficients, each coefficient 
being given by the term in square brackets in (4.29). This 
condition, then, provides an implicit relation between E and 
k. 

The matrix elements (coefficients of AA) TH. (49) 
may be evaluated to varying degrees of accuracy. In this 
approximation, since it was assumed that the atomic orbitals 
¥ centered at different lattice points overlap only slightly, 
it is consistent also to assume that the sum over j in (4.29) 
need be taken over nearest neighbours only. 

If the assumption that the atomic orbitals centered 
at different lattice points do not overlap is dropped, and if 
the overlap integrals are included in the calculation, then the 
sum over j in (4.29) may reasonably be extended to include 
neighbours other than nearest until the matrix elements converge. 


To do this equation (4.27) is replaced by 


ee ed Ge: Cee ede 2 (4.31) 
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If we now multiply equation (4.26) by Xni and integrate over 
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and again, since we sum over j we may set i = O and neglect 
the subscript to get 


A non-degenerate state now yields 
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and a degenerate state will again result in an implicit 
dispersion relation in the form of a determinant set equal 
to zero. The sum over j, now, need no longer be restricted 


to nearest neighbours. 
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Orthogonalized Atomic Orbitals 

Lowdin (1950) suggested a somewhat different procedure 
for extending the tight binding method of finding dispersion 
relations. The procedure consists of replacing the atomic 
orbitals yee by an appropriate linear combination Pq (Pr) 
of the form 


es eG): (4.36) 
J 
The coefficients W414 are chosen such that if the Aine Torn an 
orthogonal set on the subscript n (which refers to different 
atomic energy levels at a given lattice position), then the 


© Wii Torm an orithofonal set on both n and i, that is 


ni. 


fea.) Pm (r) dt = Bn 4 5° (4.37) 


Lowdin showed that 


a5 [(1¢s) 1/2 or (4.38) 


where the elements of the matrix 5S are given by the overlap 


integrals for a given level of the atomic energy E.3 thus 
g(n) = a a a ee (4.39) 
ne ae tr 1d 
; a1/2 ; canted 
The matrix (1+S) is defined by the expansion 


(4s) 1/2 = r - (1/2)8 + (3/8)8* - (5/16)89 +..., 
(4.40 
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Since the x's are atomic orbitals which are at least reasonably 
localized near a given lattice point, the expansion (4.41) 
should converge quickly to give the necessary coefficients 
ek A knowledge of the W14 permits the construction of the 
orthogonalized atomic orbitals 4 from equation (4.36), and 
these functions are then used in place of the Kew in the 
formalism of the tight binding method. Slater and Koster 
(1954) proved that the Lowdin functions 4 Mave symmetry 
properties corresponding to the atomic orbitals from which 
they were derived. Ans state, then still leads to an explicit 
dispersion relation, and in general the order of the secular 
eqaution supplying the (E,k) relation is the same whether one 
uses atomic or orthogonalized atomic orbitals. 

One requires, of course, the elements of the matrix 
(a+s) 71/2, a general procedure for finding these has been 
given by Lowdin (1951). 

We note here that in order to increase the effective- 
ness of both the tight binding and orthogonalized atomic 


orbitals procedures, on could extend the sum 


Lar (r) a iy Vin (F) (4.20) 
n 
to include terms of different energies. One could then use 
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a variational technique to determine the energy for a given 


k. 


Linear Combination of Atomic Orbitals - Interpolation Procedure 
In the two preceeding sections we have seen that when 


tecrystal orbital VE is approximated to by a linear combination 
offatemic orbitals one is led to insist that a determinant 
involving E and k be zero in order that not all the coefficients 
in the expansion for Wz(r) be zero. Equations (4.29) for the 
appropriate values of m may be rewritten as 


|m(x) ot Thee Dy | (4.42) 


k 
where 


™ ° 4U.ny J. i 
M (k) = E°6 -)e i [xn (W5-V)x, 97. (4.43) 
Slater and Koster (1954) suggested that equation 
(4.42) could be used as an interpolation aid rather than as 
a dispersion relation. In this procedure on would treat the 
quantities En and Cae defined by 


neh = { x,(W4-¥) Xn 4 dt (4.44) 


as parameters, rather than as known constants. E(k, ) must 
now be determined as accurately as possible by other methods 
at enough points k, in k-space to determine the parameters 
Be 
relation which fits exactly at the predetermined points E(k, ). 
We note here that any of the various procedures which 


yield dispersion relations, as opposed to isolated points in 


and ht ah Equation (4.42) now provides a dispersion 
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(=,k)-space, could be psed in ‘the above manner to provide 
dispersion relations passing through points E(k, ) found by 
some other method. 


Plane Wave Method 


Since it is known that We(r) can be written in the 
form 


we(#) = et Pu (s) (4.45) 


where u;(r) has the period of the lattice, perhaps the first 
procedure that comes to mind is the expansion of ug(r) ehamee! 
series of plane waves. If this is done, then Wz(r) may be 
written as 


Wel) = ) c(h) of RHEL), (4.46) 
where kK, is a lattice vector of the reciprocal lattice. One 
procedure for determining EX and the coefficients C(ktk, ) is 
given by Pincherle (1960). 

— We will require the orthogonality of the plane waves 
ems 


e over a Wigner-Seitz cell, and this is easily demonstrated 


as follows. The e*“i'", being solutions to the Helmholtz 

equation with periodic boundary conditions on the unit cell 

are clearly orthogonal to each other if they belong to differ- 

ent eigenvalues. If they belong to the same eigenvalue, the 

product of the two plane waves is a third solution to the 

Helmholtz equation, which must therefore be orthogonal to unity. 
If we expand the crystal potential V(r) in a Fourier 


series to get 


v(#) -) wk) nn (4.47) 
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and substitute this along with (4.46) into Schrodinger's 
equation we find 


i 
"2 
+ ») C(R+k.) w(k,) et (Eth, +k,)-r 
La 
= 2) o(RHE) eilktk,)-r . (4,48) 
i 
-(K+k_)°r 
If we now multiply by e n and integrate over a Wigner- 


Seitz cell we obtain 


[ fetie,[° - B] o(e+k,) +) c(R+R,) w(k,-§) = 0, 


2 (4.49) 


That is, for non-trivial solutions C(k+k, ) we required the 


vanishing of the determinant 


|™ - EI| = 0 (4.50) 
where 
Ms = wk, - k.) _ te 3 
= [k+&,/°+w(5), 1 3 (4.51) 
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It is easy to show that the standard variational 
procedure when applied to (4.46) also leads to conditions 
(4.41) and (4.42) for the determination of E(k) and C(k+k, ). 


Orthogonalized Plane Wave Method 


The potential V(r) under cmsideration in Schrodinger's 
equation is a crystal potential which is often written in 
the form 


V(z) -)v (2-R, ) (4.52) 


where v (2) may be termed a "site" potential. If ‘iv is indeed 
an atomic potential, then an expansion such as (4.46) of W in 
terms of plane waves, or any other complete set of functions, 
must converge to the energy of the innermost core electron 

if enough terms are used. (See Herring (1940)). Since these 
are not the states with which we are concerned in energy band 
considerations, steps must be taken to prevent this from 
occurring. Herring (1940) developed the orthogonalized plane 
wave method to overcome this 3 I ee a 4 9 


Briefly, the procedure consists of using as basis 


i { 3 
functions, instead of the functions gk ™, the functions 
v, (Kz) = et(etk,) +r y (3) 
oe Mag Wee ts) 
where 
=} etKFn y (2h ) (4.54) 
PRs 8 n ed 
n 


and 


[andtte tiv be fehrete Sit Bard — of vas 

ano ts tbroh of absest conte (3H.8) ot bet, a fark 

.(,3aa1)9 bas (3)2 to nobteateress\ ‘ett “ot (st.8) 
Sie 

: bodeM sve omer? 

s'lteyunibordem® nt nottseuseblenoo Tebtv (a)¥ istinetog « 

taw sette et dotdw felinedoq a set Y sobtsupe 


(92, ( 9-4) "y (= (a)v . 
+ 
part St V I. tek inetog “site” s panies sd yam (' "Vendy 
wt ¥ to (Ab.4) 28 Nove nefansaxe oe mont fettnsieg ra} ms: 
,aneliTomr) io 2982 sfeignos  taivo YNs to able one lg 
coptesls 910° se0mtegrt edt to yarene ors ot 2 


ped sone. (( ONL) Senaeae 6s@) . beau oa pans 
hisd verses ia heareonos sis ¢w dotdw at hw eotsia on 

morl eatdt travenq ot: nedist ed teum egede: a pir 
stata Ssstisnogonsse ond beaoloveb (oveL) 
¥ + LISEVELD ett enooreve "8 


atead 985 i nes to ateterop aheaget ae 4 
. | POAT, 


66 


i(k+k;.)-r _*,- 
May - fe (tk; ) or ¥y(R) adts (4.55) 
J 
The functions X.,(P-R, ) are atomic orbitals, assumed 
to be orthogonal and non-overlapping, that is one assumes 
that 


[xgFR,) xX, (r-R, ) qT = Db pea (4.56) 
aie serve tal orbitals 92, (r) are the core states one wishes 
to avoid, and the author has shown that 


ee { 
[i 2) 9g, (r) dt = 0 (4.57) 


_! 
for any values of k and s. 
One may now set 


v(B) =) 8, (ke) ¥, (BF) (4.58) 


and use a variational procedure to obtain the secular equation 
and thence an approximate dispersion relation. 

We note here that if v (2) in equation (4.52) is a 
site potential which includes the effect of the core electrons, 
then the plane wave expansion (4.46) for Ww will not, as 
contended by Pincherle (1960) converge to the wave function 
for the innermost core electron. However, one would still 
expect an expansion of Wy in terms of orthogonalized plane waves 
to be better than a plane wave expansion (Callaway (1958)). 

The plane wave expansion must rely on the potential and the 
variational procedure to make the final result orthogonal to 
the core wave functions, whereas the orthogonalized plane 


waves already satisfy this condition. Fewer terms in an 
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orthogonalized plane wave expansion should therefore be 


required to achieve a given level of convergence. 


Augumented Plane Wave Method 
Another procedure for constructing a suitable set 


of basis functions in which to expand V(r) was developed 
by Slater (1953) and Saffren and Slater (1953). With this 
method the Wigner-Seitz cell is considered in two parts, one 
a sphere surrounding the ion at a lattice point, the second 
part outside this sphere but inside the unit polyhedral cell. 
The radius ro of the sphere is arbitrary, except that the 
sphere must be contained in the unit cell. 

Outside the sphere the potential is taken to be 
constant, and inside a spherically symmetric approximation 
is used. The basis functions are constructed in the 
following manner. 

_ _ Outgide the sphere every basis function is given 
by et Ws thas plane wave is written in the form 
c= 4 
olKF =) a (aett) P, (cos 6) J, (er), (4.59) 
y} 


where the P ,(cos@) are Legendre polynomials and the j (xr) 
are spherical Bessel functions. 9 is the angle between k 
and r. 
Inside the sphere the different basis functions are 
given by u(r,@,E) for values of E yet to be determined. u(r,@,E) 


is given by 


u(r,@,E) -) a, u, (r,E) P (cosé ) (4.60) 
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where u ,(r,E) is the solution for a particular E, finite at 


the origin, of Schrodinger's equation with magnetic quantum 
number zero. 


To make the basis functions continuous, the coef- 
ficients a, are chosen to be 


\ 


a, = i°%(241) 3, (kr,)Ay, (ro). (4.61) 


The derivatives, however, of the basis functions across the 
surface of the sphere are discontinuous. 
If we now set 


Yz(r,6,E) = u(r,9E), O ae o (62) 


then the set of basis functions is given by Y,(r.6,E,). The 
set. of values (E, ) are chosen to be the stationary states 
With respect to E of 
% ; : , ‘ 
<2 = [Ss (r,0,E) H YE (e,0,8) dt (4563) 
It is shown by the authors that the stationary states Ey oP i 
< E> are obtained when E = Ey and that 


E, = «2 + (4or,-/v) ) (2044) (xr) < Inu ,(roE,), 

: (4.64) 
where here v is the volume of a unit cell minus the volume 
to the sphere of radius r,. 

Equation (4.64) must now be solved for E,, i = 1,2,3,.. 
and ths basis functions Yz(r,9,E,) constructed and used in gay, 


a variational procedure. 
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CHAPTER V 


NEW METHODS FOR CALCULATING ELECTRONIC ENERGY 


BANDS AND FERMI SURFACES 


In this chapter we describe first a fairly standard 
matrix method of finding approximations to the solutions of 
a differential equation, and we consider in some detail the 
usual variational procedure. We then develop the "restricted 
cellular basis function" approach, which is essentially an 
extension of the method of "solid harmonics" suggested by 
Betts (1959). A second new method based on a cmsistency 
consideration is then described, and finally an interpolation 
scheme based on the two new procedures is considered. 


Equivalent Matrix Equation 
Variations of the procedure developed below are not 


uncommon; see, for example, Bickley and McNamee (1960). 
Let us consider the equation 


Hy = © (Sui) 


Where H is a linear differential operator and y is to satisfy 
the same homogeneous boundary conditions as does the given 
function ®. If we have a sufficiently complete set of basis 
functions (f,) satisfying the boundary conditions of (5.1), 
then we may write 


Vv => ee . (5.2) 
ro a 
ds 
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To look for an approximate solution to equation (5.1), let 
us assume that we can write 


CO 
= af. .3) 
v ) aan (5.3) 
J=1 
and 
: oe) 
@ = core. (5.4) 
) egf 5.4) 
j= 
The coefficients a Bee (5.3) are to be found, and the 
coefficients C , in (5.4) are defined by 
Ne = 6 (f65) 


Where, of course, the i'th component of ¢c is c,, and the 


matrix N and the vector © are given by 
N.. = /f£.f.dt (5.6) 


and 
(6), _ [eae (5.7) 


Substituting the expansions (5.3) and (5.4) into the differential 
equation (5.1) yields 


) asf, =) ety. (5.8) 
J J 


If we now multiply equation (5.8) by f, and integrate, we get 
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! 
the i th row of the matrix equation 


Ha = Ne (5.9) 


where H is now a matrix, the elements of which are given 
by 


Ay = Jeet dt (5.10) 


and the components of the vector a are the as. Since we can 
multiply (5.8) by any f,, i=1,2,...n, the entire matrix 
equation (5.9) holds. 

An approximate solution to (5.1), then, is given by 


(5¢3); where the coefficients a, are found from 


i 1 


a=H “Ne = H-6 (5.11) 


This procedure is easily modified to treat the 
eigenvalue equations we have been considering. If we put 


6 = Hy, (5.12) 
then it is clear that equation (5.9) is replaced by 
Ha = ENa, (5.13) 
The eigenvalues and eigenvectors of (5.13), then, provide 


approximations to the eigensolutions of the differential 
equation Hy = Ey. 
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Variational Procedure 


Consider the ditferential equation 
Hy = Ev (5214) 


plus homegeneous boundary conditions; if we multiply by v 
and integrate over the appropriate volume (S14) leads to 


: __feavar : (5.15) 


fon 


It is easily demonstrated as follows that if H is Hermitian, 
then those functions Ww. which produce stationary states En 
of E (as defined by (5.15)) with respect to variations in v 
are identical with the eigensolutions (EV) of equation 
(5.14) 

To do this let us expand yw in terms of a sufficiently 
complete set of functions (f,), and then vary W by varying 
the coefficients in the expansion. Let us say, then, that v 
is given by 

00 
Yom ae (5416) 
Leak 
Interchanging the order of integration and summations, 
substitution of (5.16) into (5.15) leads to 
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where we again define 


Hy , - [ eyue jac (5.10) 


and 


Nay - f ey2,a0 (5.6) 


Equation (5.17) may be written then in the form 


(5.18) 


Where H and N are matrices with elements Hy and Nay 
respectively. To find the stationary states of E we 
differentiate with respect to a = eos sekeew ne She2eo 


that the results vanish. 


\ 


oe 
_ a:(Ha dan, +) ame, | = 0 
a : ; (5.19) 
That is 
7 [aTana, i oe [ (x7), or (5.20) 


and since this must be true for all k we find 


(8L.2) 
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re ee a: (H ee 
(H'+H)a = a (oan) a. (6.21) 
‘Now, N is clearly symmetric, and go is the matrix H if the 
differential operator is Hermitian; withthe aid of equation 
(5.18) then we may rewrite (5.21) in the form 


Ha = ENA. (5,20) 


But according to the previous section this matrix equation, 
where H and N are now infinite square matrices anda is an 
infinite column vector, is equivalent to the differential 
equation Hy = EY plus boundary conditions. The stationary 
states of (5.15), then, are the eigensolutions of (5.14) and 
alee of (5.22). 

We note here that if the differential operator is 
not Hermitian, then (5.22) must be replaced by 


(1/2) (H+) = ENa, (5.23) 


and we have no guarantee that the eigensolutions of this 
equation bear any relation to those of HW = EY. If the 
differential operator is Hermitian, then the variational 
procedure leads to exactly the same result as does the 
equivalent matrix approach. However one important advantage, 
outlined below, has been gained. 

To use the variational procedure as a means of finding 
approximate eigensolutions to Hy = EY, we merely use a finite 
Subset, say the first n, of the set of functions (f,), instead 
of using the whole set. But this is equivalent to insisting 
that 


a 0 ee eo oo (5.24) 
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Now, from (5.18) it is clear that E(a) = E(ca) for any 
nonzero constant c, and therefore E(a) is finite as long as 
a # O, since N is positive definite. Therefore, if there 
exists a smallest (algebraically) value of E, say E, ; which 
is a stationary state of (5.15) or (5.18), then E, must be 
a minimum value of E(a). Any restriction, then, such as 
(5.24) on the a. can only increase the smallest eigenvalue. 
Thus, the solution to 


Ha = ENa (5.25) 


where H and N are now nxn matrices and a ig an n-vector 


provide approximate solutions to 
Hy = Ev (5.14) 


and in addition, an upper or lower bound to the smallest or 
largest eigenvalue respectively of (5.14) if such exists. 

Let us now consider the solution of equation (5.25). 
Normally one finds numerical solutions to matrix equations 
of the form 


where A is a matrix, y a vector and A a scalar. There are 

a variety of ways in which (5.25) may be reduced to this form. 
We note first that if the set of basis functions Cre, were 
orthonormal, then N would be equal to I, and (5.25) would 
already be in the form (5.26). If this is not the case, one 
obvious method would be to find the inverse of the matrix 

N, and then look for the solutions of 
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N-Ha = Ea. (5.27) 


This procedure, however, should not be used since it involves 
considerably more labour than the well knownmethod described 
below. 

Since N is positive definite it is always possible 


to write 
N= mM (5.28) 


where M is uniquely defined by this equation if we insist 
that M be lower triangular, that is, M have only zeroes above 
the principal diagonal. M is relatively easy to find 
numerically, and equation (5.25) becomes 


Ha = EMM‘a, (5.29) 
or multiplying on the left by M+ 
mtHa = EM a. (5.30) 


Let us now define the vector 


in which case 
a= (m*)te , (5.32) 


and (5.30) now becomes 


(78.2) 
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| mw n(n)? 5 = Ec (5.33) 


which is in the required form. The eigenvalues of (5.25) can 
now be found from (5.33) and if necessary one can also find 
c and thus a. 

We note here that mt is Yéelatively easy to find, 
since Mis triangular. The method described above is equivalent 
to using the Schmidt orthogonalization procedure to replace 
the original set of basis functions by a new orthonormal set. 
Assuming orthonormality usually results ina simplification 
of analysis. As a practical procedure, however, it is usually 
more convenient to use non-orthonormal basis functions, and 
then resort to (5.33), rather than (5.25) with N = I. 


Restricted Cellular Basis Functions 
The standard cellular methods described in Chapter 
IV, for solving HW = EW where H contains a periodic potential 
all have one feature in common; they all require that the 
potential V(r) be approximated by a spherically symmetric 
potential inside some sphere. The sphere may be an approxi- 
mation to th entire unit cell, as in the Wigner - Seitz 
procedure, or it may be a sphere inside the unit cell, as 
in the "lattice harmonic" method of Slater. It has been 
demonstrated by Betts (1959) that this approximation can in 
some cases lead to serious errors; the more the unit cell 
deviates from a sphere and the fewer the group operations 
under which the wave function is invariant the more likely is 
it that the spherical approximation will lead to inaccuracies. 
The cellular method described below overcomes these difficulties. 
Now, given the equation Hy = Ey we are looking 
primarily for E(k) where ¥ = e "Tul (P). Basically the method 


| ‘i 


(€€.2) 38 1 fhe nyete | 


ms (2S.¢) to esslavnegie anit saris? bert 
baft os fe aga SHO yrssasoen 2 bris (8€-€) 


<bntt ot Yess ylovitéler ar 8 ‘ut oni oben 9 
jnsisviups e@f svods Bsdbtogsb Fortean eiT Laeivedsiasy & 
sosiqer oF srubsoorg fo Mest isnogon tio: thimiioe oat wy 

.te% [sortenodinro wen B.y0 ahotyonu alesd to tea Lf, 
aotisolttftumta s 7 at inest vitsver vd tLewroriods 29. 
(ifsray eb tL ,tevewor! (saebecogg issttosnq we ah 4 é. 


bis ene tion atesd Lemenodd10-nex samt 


l= dtiw (as.@): tents tediet <(aneeh 
: 08, 2S 


satel 


~ POLS S029 ELE Se 
isvgeto ot bedireesb ehonten telvifec bt waa 
fatiasi0qg olholusg & antsines H se tertw a = « ot “ 
sit sady ettupes Ifs yeds s:nomndo @ ; rs 
ol td ominrye viisetvsiza 8 Yd betsmixorwggs od (a)¥ £ £ ined 
-~txotgge ns sd vem eratge eit .sierige emon sbtent £ tere 
stis® - rempiW ed? at as .ifen tine erlgae ot Daas 
as .[ieo stay sdt ‘Sblaot eretqe so ed yao tf 40 | 
esd ash 2 .tetel® to hottem “otsomend sotdd. 
ut iiss noitembxorggs ates dents (eet) hacabeotetie 
[Leo tihw eng eon edt LetorTTs score ; 
enotvs eqs gory ent ‘rows? ons Bae | 
et yell etom ent, tasia 
-aaboetiucosnt or best Citw 
.eetclvoltitb Sead semosasve © 


ep: 
| ois 


78 


consists of expanding ¥ in terms of a set of basis functions 
flr); and then using the variational procedure described in 
the previous section to determine the energy for a particular 
value of k. Any periodic potential may be used without 
resorting to a spherically symmetric approximation, and the 
boundary conditions, as supplied in part by the Bloch form 

of ¥, are satisfied exactly on the surface of the Wigner- 
eeitz cell. 

Having decided to look for E(k), one first requires 
that subgroup G({k) of the point group G of the lattice, the 
operations of which either leave k invariant or change k into 
an equivalent reciprocal lattice vector. Since the irreducible 
representations of G(k) determine the transformation properties 
of VitR =e Uny these irreducible representations must now 
be found in order that the possible symmetries of ao may be 
determined. Here the range of n is the number of such 
irreducible representations, and we note that except for 
special values of k on lines or planes of symmetry in the 
first Brillouin zone, G{k) will be the identity operator only. 
In this case the number 1 is the only irreducible repre sentation 
of G(k), and we find only the uninteresting identity 


We, (Ir) = Wg, (7); (5.34) 
“that is, W has no special symmetry properties. 

Having determined the possible symmetries of ¥ 
(assuming that k is ona’ line or plane of symmetry) on must 
now decide to which of the irreducible representations wy is to 
correspond. In general, the greater the symmetry of vw, the 
smaller the corresponding eigenvalue. Once the irreducible 
representation has been chosen, the symmetry conditions wich 


wv must satisfy may be tabulated. 
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Given the symmetry conditions of WwW it is now 
possible to determine more stringent boundary cmditions 


which wW must satisfy than those supplied by the Bloch 


conditions 
- a ae 
VAT a) (3535) 
and 
ay(7,) _ _, May(z,) 
= ae (5.36) 
- 2 


Here an and t, are corresponding points on opposite faces 


of the Wigner-Seitz cell; that is 


where R is a lattice vector. 
For example, if k = 0 equations (5.35) and (5.36) 


become 
W(r,) = W{r,) (5.38) 
and 
av(r,) _ _ awlry) (5.39) 
on. an, 


that is, Ww must be periodic with the period of the lattice, 
which is exactly the expected result. Now, if the group G 


of the space lattice contains for every boundary point ry 
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an operation B, with corresponding rotation matrix S such 
that 


Sr, = fo, (5.40) 
then 
g wr) _ av(Sr,) _ av(rs) | (5.41) 
a a(Sn,) aN, 


If, in addition we are considering the identity representation 
or G, then 


pC Sia a, (5.42) 
on, an, 
therefore 
oa — 
ay ee! (5.43) 
on, on, 


Equations (5.39) and (5.43) together imply that 


ay(r) = 0 (5.44) 


én 


for any vector r on a face of the Wigner-Seitz cell. 
On the other hand, cmsider a representation of the 
group G(k) in which 


et (5.45) 
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for all elements of the group. In this case, for those 
points ry and Ps for which &. = +1, condition (5.44) applies; 


for those faces of the unit cell for which a. = -1 if 
Sry = Ue then 

ew(r) = v¥(Sr,) = ¥(r,), (5.46) 
and 

gv(r,) = (-1)¥(r,), (5.47) 
that is 


\ 
} 


Using this result together with equation (5.38), we find 
that on these faces of the Wigner-Seitz cell 


w(r) = 0. (5.49) 


For values of KAO; and for irreducible representations 
of degree more than unity, more complicated, though less 
stringent, boundary conditons may apply. In extreme cases, 
for example if G(k) consists only of the identity operator, 
then the Bloch conditions themselves may be the most stringent 
boundary conditions available. 

Given the symmetry and boundary conditionsthat Ww must 
satisfy, the basis functions fey may now be constructed. 

In principle the BA) may be taken as linear combinations of 
any complete set of functions, the linear combinations being 
chosen such that each oy 2) satisfies both the symmetry 
conditions and the boundary conditions on Ww; clearly then, vy 
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which is a linear combination of the f'ts, will also satisfy 
these conditions exactly. In practice ome would use either 
powers of the coordinates, or the eigensolutions of a 
Hermitian differential operator to ensure completeness of 
the original set of functions. In Chapters VI and VII we have 
used only the powers of x,y, and z, and trigonometric functions 
of the appropriate periodicity, that is, the plane wave sol- 
utions of the Helmhotz equation. 

Construction of the basis functions AP) is not 
an easy task. For polynomials the procedure consists 
essentially of writing down the general polynomial 
j_k n 


a a ele Z 


P(xsy.2) =a ta xy ae 


arr = ee . 
S00: ~ 100 ijk 


for a sufficiently large n; the application of the symmetry 
and boundary conditions to Pa will result in a large number, 
M, of linear homegeneous equations in the N unknowns Bs a" 
Since an infinite number of linearly independent polynomials 
must exist which satisfy a finite number of consistent 
conditions, we are guaranteed non-trivial solutions. N, then 
must eventually increase faster with n than does the number 
of linearly independent equations, resulting in more unkowns 
than equations. It has been found, in fact, for sufficiently 
symmetric lattices (in particular, plane hexagonal and face 
centered cubic) that the M equations will not be linearly 
independent, and thus solutions can be found before N>M. 

The computational problem is by no means trivial. 
For example, in looking for a polynomial P satisfying only 
the conditions that P have the same value, and @P/an have 
the same magnitude with opposite sign at corresponding points 
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on opposite edges of a plane hexagon, it was found that the 
only such polynomial of degree four or less was a constant. 


Setting 

P = agox°ta x-y+...ta yta, (5.51) 

51 O1 OO 

resulted in thirty three homogeneous equations in the twenty 
eight unknowns 25 5k" These equations have three linearly 
independent solutions which may be chosen to be of degree 
zero, five, and six respectively. The use of automatic digital 
computers in such circumstances as these is almost mandatory. 

The difficulty involved in finding the basis functions 
ace) is not without its compensations. Firstly the f AD), 
once found, can be used with a variety of potentials and 
therefore for different materials having the same structure 
or with different trial potentials for the same material. And 
secondly it is not unreasonable to assume that since the 
boundary and symmetry conditions place such stringent restrictions 
on the f,(r), that relatively few such basis functions should 
be needed in the variational procedure to provide a given 
level of convengence in the energy. That this is actually 
the situation at least in some special cases is demonstrated 
in the following chapters. 

Let us now consider the use of trigonometric functions 
of the appropriate periodicity to construct the basis functions 
EAD) The trigonometric functions are solutions to the 
Helmholtz equation, and therefore the empty lattice problem. 
According to Chapter II, then, one can find symmetrized 
linear combinations of the degenerate states which transform 
according to the various irreducible representations of 
G(k). But since the boundary conditions which ¥ must satisfy 
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do not depend on the particular form of V(#), only on its 
periodicity, then these symmetrized linear combinations of 
the trigonometric functions must already satisfy the boundary 
conditions of Hy = EV with VO, as well as the empty lattice 
problem. With trigonometric functions, then, the determination 
of the neo) is considerably simplified; once the symmetry 
conditions appropriate to the choice of k and the irreducible 
representation have been applied, the boundary conditions 
must automatically be satisfied. 

We note here that if k is not chosen on a line or 
plane of symmetry of the first Brillouin zone, then the 
method described in this section reduces merely to an application 
of the variational procedure. In particular, if the Pm) are 
constructed from trigonometric functions then the procedure 


becomes equivalent to the Plane Wave Method. 


Consistency Determinant Approach 
The method described below, which is similar to 


Hill's procedure for treating the Mathieu Equation ( Hill 
(1886), McLachlan (1947)) is based on a cmsistency requirement 
on ug(P) where 


ve(@) = e'* Pu, (zr). (5.52) 
We consider 
[-V°sv(2)]v = By faoea) 


Where V(r) is periodic; if we substitute (5.52) into equation 
(5.53), we find that ug (Pr) must satisfy the equation 
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ee oe ee 
Vout2ik-yu+(E-k-k-V)u =0. (5.54) 


Let us now separate ug (7) into its real and imaginary parts, 


and write 


u(r) = e(r)tiy(r) ; (5.55) 


substituting (5.55) into equation (5.54) yields the two real 


equations 
= Fis ee eee 
Va -2k-VY+(E-kek-VS d= 0 
and 
2 i. 


Vo v+2k-va+(E-keK-V)y = 0. (5.57) 


We now expand both @(r) and y(r) in terms of a set of basis 


functions PAP): 


ate) i=) dat (2) (5.58) 
and 
y (2) =) cage. (5.59) 
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Substituting the expansions for @(f) and y(r) into equation 


(5.56) gives 


avr, 5-2) cskeve ot B-R-k)) oe Bee (5.60) 

me J sth )) 238 s-) at, ' (5.60) 
J 3 

=. us now define the matrices L,N,A(k), and V by the following 

equations 


by = [Pye (aya: (5.61) 

1368) = [t,@)R-ve ,(#)ar (5.62) 
ee core (5.63) 

Vagve is (r)V(r) )£,{ r)dt (5.6 


Where each of the integrations is over the Wigner-Seitz 
cell. If we now multiply equation (5.60) by f (r) and 
integrate over the unit cell we find, using me definitions 
(5.61) to (5.64), the 1 th row of the matrix equation 


(L+(E-k-k)N-V)a-2A 


Pa 
yl 
ed 
QI 
H] 
a 
\e 


(5.65) 


Where the coefficients as and Cs comprise the vectors a and 
Cc respectively. Since we could have multiplied (5.60) by any 
of the basis functions Eee r it is clear that the entire 
matrix equation holds. 

Similarly, if we substitute the expansion (5.58) and 
(5.59) into equation (5.57), and then multiply by f(r) and 
integrate, it is easy to show that 
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(I+(E-k-k)N-V)c+2a(k)a = 0. (5.66) 


Clearly, then if we want non-trivial solutions for the vectors 
a and c, then equations (5.65) and (5.66) imply that the 
determinant of the matrix C(E,k) must vanish, where 


It+(E-k-k)N-V ~2A(k) 
C(E,k) = (5.67) 
+2A(k) L+(B-k-k)N-V 
If we use M basis functions fAE Ts then C{fE,k) = 0 provides 


an approximate (E,k) relationship, where C is a (2M)x(2M) 
matrix. 
The labour involved in evaluating determinants of 
large order is not insignificant, and therefore the MxM 
formulation derived below is considerably more convenient. 
Let A and B be arbitrary MxM matricies: we consider 


now the problem of reducing the 2MxeM matrix 


A -B 
(5.68) 
B A 
to the block diagonal form 
C O 
(5.69) 
O D 


where C and D are also MxM. We treat the problem much as if 
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A and B were scalars. 
Let a and b be scalar multiples of the MxM unit 


matrix I, and assume that L is an M x M matrix such that 


A -B a a 


B A b b 


The matrix equation (5.70) may be rewritten as the pair of 
matrix equations 


(A-L)a-Bb = 0 Gras 
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Solving both (5.71) and (5.72) for {A+L) we find respectively 


A-L. = Bba (5.73) 
and 
_o-L fo om) 
het Sean 3 (5.74) 
therefore 
ae) ae (5.75) 
Let 


bat = Cy, (5.76) 


ay a ra, af — 8 
tinw MeM sit to eotgnanua ‘aeispe ad _ sod 
tend? dow xitten M “7 nis: BL J sed. mm 


io “ieq edt 6 Io7 sesh od yen (OF. 2) mish 2 eb 
Sno 


4 


0 = da+2{ I-A) 


‘it 
a0 = OL. A)vee 


s 


vitosqser batt ow (268) aot (ot, =) 


Ppt x oak 


< dea- = J.A 


at ae : 


where the matrix c is also a scalar multiple of the MxM 


unit matrix I, since both a and b are. From (5.76) 


but from (5.75) and (5.76) 


therefore 


and so 


And since c is a multiple of 


jy = 
Let us set 
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where we use either the top or the bottom row of signs. 
Equation (5.70), together with equations (5.82) to (5.84) 
now yields 


A -B a A -B iN 
B A db B A Fil 
= A+iB 
B+iA 
= si 
(A+iB). 
+il (5.85) 
Therefore 
A ~B T iL a Zz A+iB O 
B A -il me -il niet @) A-iB 
(5.86) 
and since 
a ae tie fo at 
- 1/2 (5.87) 
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we find that 
T iT 


(1/2) ‘ 


Therefore 
A -B A+4B 
A = 7 
B A 
= |a+iB| 
And since 
(A-iB) = (A+iB)* 


implies that 
[A-1B] = |A+iB] * 
where * means complex conjugate, then 
A = |det(atis)| *. 
that setting 


We find, then, 


c{E,k)| = ¢ 


A+iB 
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(5.92) 


(5.93) 
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where C(E,k) is given by equation (5.67), is equivalent to 
setting 


[D(E,z)| = 0 (5.94) 
where 
D(E,k) = L+(E-k-K)N-V+21A(k). (5.95) 


Since, according to Blochts Theorem, ug(r) is 
periodic with the period of the lattice, then the Pi\e) used 
to expand the real and immaginary parts of u must be periodic 
or at least satisfy periodic boundary conditions on the 
surface of the Wigner-Seitz cell. In this case it is clear 
that the matricies L,N and V are symmetric and the matrix 
A(k) amtisymmetric, and therefore D{E,k) is Hermitian. 

Given that D(E,k) is Hermitian, there exists a matrix S such 


that g~lp(E,&)s is both diagonal and real, and since 


Is“p(z,&)s| = |s7||D(z,e} | s] = |D(z,.%)| (5.96) 


for any E and k, it is clear that the values of E, given Ky 
which satisfy (5.94) are, as expected, real. Futhermore, 

in the special cases treated in the following chapters it has 
been found that the elements of D/E,k) are not complex, but 
are either real or imaginary in such a fashion as to permit 
the multiplication of alternate rows and columns by 1 to 
obtain a determinant consisting entirely of real elements. 
This greatly facilitiates the evaluation of specific points 
(E,k) from the implicit approximate dispersion relation 
ID(zE,z)| = 0. 
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Interpolation Procedure 

The methods described in the two proceeding sections may 
be combinedin the manner outlined below. 

The method of restricted cellular basis functions 
should provide, once the potential and the basis functions 
themselves are available, accurate values of E(k) for only 
a modest outlay of labour. One must evaluate the integrals 


2 ae 
ane [ay f dt (5.97) 
Va; = [eyve yr (5.98) 
and 
Nyy = ff jr, (5.99) 


and then the lower eigenvalues of the matrix equation 
Hy = (-L4V)y = EN (5.100) 


must be found. Standard library programs are usually available 
to solve equation (5.100) on a digital computer, as are 
integration routines requiring only the writing of an integrend 
subroutine, for the evaluation of the matrices L,V, and N. 

On the other hand the initial determination of the basis 
functions fr, 
expected to be convenient for those values of k on lines or 
planes of symmetry in the first Brillouin zone; the more 
symmetric the point k, the more severe the restrictions on 
the basis functions fis and therefore the fewer the number of 
basis functions required in the variational procedure. 
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The determinental method described in the preceeding 
section, on the other hand, provides one with an approximation 
to the entire dispersion relation. In this case, the basis 
functions being much less stringently restricted are much 
easier to find. But for exactly this reason one would expect 
to require a correspondingly larger matrix to provide a 
given level of convengence. 

It seems natural, then, to follow the example of 
Slater and Koster (1954) and combine the methods described in 
the two proceeding sections. The consistency determinant 
provides an implicit dispersion relation involving, besides 


E and k. the quantities 
0 
‘[tyae jr (5.101) 


where Q is any one of 1,V°>V; or V. These latter quantities 
may be treated as parametrs; they may be evaluated exactly, 
or by using @ least squares technique, to ensure that the 
interpolative dispersion relation fits at the symmetry points 
(E,k) determined by the method of restricted cellular basis 


functions. 
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CHAPTER VI 


APPLICATION OF NEW METHODS 


In this chapter the three new methods developed in 
Chapter V are applied to various models. The applications 
are intended to illustrate both the use of the methods and 
the rates of convengence in some special cases. 


The Empty Lattice 
The (E,k) relationship, as V(x) approaches zero, 


becomes 
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where V(x+2n) = V(x). In this section the approximations 
TO (6.1) are found by using the consistency determinant with 
various types of basis functions. 

If we use as basis functions 


f, = 1N 2 
Pon = gin nTx 

= Tx, n> ; 
Pont] cos nrx; nei (6.2) 


then the following results are easily obtained 
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bs id= : 
-[ (m-1)7/2]°, m oda. (6.3) 
Nom = Onm (6.4) 
Son, out = er 
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Ted RO (6.6) 


Here €;, - [ty lat ,/ax)ax. 

Using only the first basis function, the matrix 
aa N(E-k“) -V+2iké ] becomes (E-k") and we find exactly 
the lowest branch of the dispersion relation in the first 
Brillouin. zone, -1/2<k.< 7/2. 

Using the first two functions, we find 


(tiien Gia es (6.7) 


mak 18 2. < Ke Which is the correct lowest branch of the 
dispersion relation, and E = Koa = which bears no relation 
to the correct relation. 


If we use the first three basis functions, we find 
2 a 2 
(E-k°) (E{k-1))(E-(kt1)“) = 0; (6.8) 


that is, the lowest three branches of the (E,k) relation in 


the first Brillouin zone are obtained exactly. 
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In general, the consistency determinant consists of 


(E-k~) in the (1,1) position, followed by n 2x2 blocks along 
the main diagonal, each of the form 


“) 


B-ko-m@1* —OmTrik 


e Wale; een 


om ik B-ko mer 


where the ‘first (2n+1) functions are used; and if the (en+2) | th 
function is used also, the last 2xe block is followed by 
E-k°-(n+1)1° in the (2n+2,2n+2) position. 

Clearly, then, using the first ntl basis functions 
provides exactly the first n+l curves in the first Brillouin 
zone, and including the (2n+2)'th basis function given only 
an extraneous root. That exact answers are obtained is not 
surprising; the basis functions used are the eigensolutions 
tothe problem. That the basis functions need be added in 
pairs is not obvious, and seems to be related to the symmetry 
of E(k), that is, the symmetry of V(x). 

Let us now consider the use of polynomial basis 


functions with the consistency determinant. If we choose 


then setting the determinant equal to zero yields 
Siege Sal 
(B-k~)* = 0 (6.10) 
when n basis functions are used. The difficulty, of course, 


is that the f'ts do not satisfy the periodic boundary conditions 
which u(x) must satisfy. If the more sensible basis functions 
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fs nD 
f, = a 
P. = xt -2x° (6,14) 


which satisfy the boundary conditions 
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are used then the following results are easily obtained. 
LtN(E-k°)-V+2iké becomes 
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2(E-k*) O - 14 (p12) 
8 + 16 2 TO ac 
Wp 42 128, _ 256 + 214i, 12 
- [5 (E-k") - To5"4 105 315¢E-**) 
(6.13) 


Expression (6.13), using only one basis function yields 


E =° as expected. Using two basis functions we find, besides 


E =, i, 


A =. 80/2 4k". (6.14) 


This is not an improvement, and would seem to correspond to 


the extraneous curve E= “rane for the trigonometric basis 
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functions when n= 2. Forn = 3, (6.13) yields E = k- and 
2(E-K°)° - 41(E-K2) + 210 - BoK2 = 0 (6.15) 


Which gives 


E = k© + 10.25 + .25(1 + 640K")2/2. (6.16) 


Figure 4, may be used to illustrate the situation. Using the 
minus sign gives the middle solid curve (-1/2< k < 7/2), and 
the plus sign gives the upper solid curve. These two curves 
form approximations to the parabolic segments E = (x-1)* 

and E = (ktr)*, -1/2 <k < 7/2. The following table illus- 


trates the accuracy of the approximation. 
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- “i ESxact al ce 
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0.5 Tios 6.98 13.67 13.20 
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1.5 3.01 2.69 21.99 21.54 
1/2 2.78 2.47 22.65 Se.ee 
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the intersecting of the two parabolas E = (ktr)°, akin to the well 
known difficulty in finding a polynomial approximation for 
|x| , and is not likely to cause trouble in a lattice with 


v0. 


The Kronig-Penny Model 


We now consider the solution to 
(-V4V)y = BY (6.17) 
where 
V(x) = 2Vd(x-en). (6.18) 


The exact dispersion relation for the problem is easily found 
to be 


cos 2k = cos WE + 2V sinwE (6.19) 
NE 
(see, for example, Jones (1960), 
Let us consider first the method of restricted basis 
functions. For the Kronig-Penny problem the point group G of 
the lattice consists only of two elements G_ and G. where 


1 2 
x = x and Gx = -xX. The two irreducible representations of 


te group are the identity representation and the alternating 
representation. Both representations apply to the only two 
points of symmetry k = O and k= 7, and so we have in total, 
four cases. The boundary conditions which the basis functions 


must satisfy are easily found, and are as follows: 
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Identity representation, k = 0 p(-x) = p(x) 
pt(+1) = 0 
(6.20) 
Alternating representation, k = 0 q(-x) = -q(x) 
q(+1) = 0 
P6ve14 
Identity representation, k = 7/2 r(-x) = r(x) 
eel) =O 
(6.22) 
Alternating representation, k = 7/2 s(-x) = -s(x) 
at(+t) ae 
(6.23) 


The basis functions for the four cases are easily found: 


p(x) = nx-"-*_(n-1)x=" (6.24) 
a, (x) a xen t_yen-3 (6 25) 
oe9, = xen-2_,en-4 (6.26) 


and s(x) (2n+1 


Let us consider the alternating representation for k = 0; 
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the other cases are treated in an almost identical fashion. 


We need first the matrices 


i 
Ls = [ay (x)a"\(x)ax (6.28) (6.28) 
-l 
al: 
Veg =f ay(x)VOx)a, (x)ax (6.29) 
-1 
i 
Nis =f ay (x)a,(x)ax. (6.30) 
-1 
The matrix V is identically zero, since q.(0) = 0, and L and 


N are easily evaluated exactly; L and N are symmetric, and 
117 -24/15, and Ni 16/105. Using 


seccessively the first 1,2,3 and 4 basis functions we find the 


typical elements are L 


eigenvalues of Ha = ENa where H = -L+V = -L. The following 
table results 


Taple ° IT 

n EL Bs 
1 10. 5OG0G m « leone t= Sees 
2 9.87539 50.12461 
9.86962 39.99799 
4 9.86961 39.48929 


ie.e) 


(exact) 9.86960 39.47842 
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That V is identically zero is not surprising; equation (6.19) 


implies that every second root E for k = 0 is given by 

ANE = nr (6.13) 
or 

E = (nr)* (6.32) 


Then E, = nm and E, = ba, regardless of the strength of the 
delta function potential. The convengence is quite rapid, 
and the alternate roots (which do depend on V) ror k = U,. and 
the roots for k = 7/2 are easily found in the same manner. 
Again, Figure 4 is a good illustration of the situation. 

Let us now apply the consistency determinant to the 
Kronig-Penny model, with V = 1/2, i.e. V(x) - ) 8(x-2n). We 
consider first the trigonometric basis functions 


f, = 12 
Pon = att ny & 
Pont] = cos nr x,n>l. (has) 


The matrices L, N and € have already been used for the empty 
lattice and are given by equations (6.3) to (6.5) inclusive, 
only the matrix V is now required. V is given by 
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and so Vii = 1/2 
Vy 003 = Yo35,1 1Ne2 
oa ae 
aJjsed 
Va; = O otherwise. (6.34) 


For n = 1 we find 


a en ae (6.35) 


For n = 2 equation (6.35) holds, and we also find the extran- 


eous curve 
R=aw +k. (6.36) 
With the first three basis functions we find 


O(a) | (Bako sa) (Be ge 1) Sk] (Gk) —_oe 


(6.37) 


At k = O the three roots of (6.37) are .45200, a. and 
10.91760, whereas the correct values are .42695, and 
10.84276. At k = 7/2 the three roots are /4 = 2.46740, 
3.44145, and 22.73256, while the corresponding exact solutions 
are 1/4 = 2.46740, 3.37374 and 22.20657. The solutions to 
(6.37) are found in tables IV, V, and VI. 
Let us now use the consistency determinant with poly- 
nomial basis functions satisfying the periodic boundary conditions 


(6.12). One such set of basis functions is 
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ts aE 
r Z xooti_,en-1 
en 
ee ees eo (6.38) 
‘onit.= 2 er eck Pe 


The matrices L, N, V and € are found to be: 


Pod lon : i ete —% 


L - -8 (8mn -2m-2n-1 
en,2m =s- (@m+2n+1) (2m+en-1 ) (2m+en-3 


-64 (mt1)(n+1 


eet eal ~ (om+en+3 ) (om+ent1 ) (2m+en-1 (6.39) 
Nj, =2 
di geg! 
-2(4n+3 


Ny ont ~ n(ent1) (ent+3 
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(6.40) 
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11. = S14 = 


é 2 -€ = -32 ntl 
ent] «em en,emt+l (emtent3 ) (2m+2n41 ) (2m+en-1 


bs4 = 0 otherwise (6.41) 
and 
Vi; =O tr a ee 
Vi, =1 (6.42) 


Using only the first polynomial basis function, we find, as 
we did for the trigonometric basis, 


Se ee, (6.43) 


‘ 


and using the first two basis functions gives, as well as 
(6.43), the extraneous root 


ec Ei iat (6.44) 


which corresponds to (6.36) for the trigonometric approximation 
to the Kronig-Penny model, and to the extraneous curves when 

n = 2 for the empty lattice. If the first three basis functions 
are used, one obtains 


5 (E-k°-1/2) (2(E-k*) -21) (107 (E-k~) -384) 


~-34.3(E-k°)*(2(E-k°) -21) -15360k"(E-k#1/2) = 0. 
(6.45) 
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At k = O the three roots are .45720, 10.5, and 10.93605 

as opposed to the exact solutions .42695, a = 9.86960 and 
10.84276; at k = 7/2 the roots of (6.45) are 2.60478, 
3.57879 and 23.11187, while the true roots are 2.46740, 
3.37374 and 22.20657. The solutions to (6.45) appear in 
tables IV, V and VI. 

Using the first four basis functions produces only 
an insignigicant improvement in the three lowest (E,k) curves; 
the results for n = 4 agree with those for n = 3 to almost 
three figures. This corresponds to the empty lattice situation 
in which the symmetry of the problem appears to demand that 
the basis functions be added in pairs to improve the approx- 
imation. We would expect, then to get an improvement for 
n= 5, and this is indeed the case, as illustrated in the 
following table. 


Table III 

n E, (0) E,(0) E(0) 

if 5 7 ; 

3 ~45720 10.5 10.93605 
5 44625 9.87539 10.89079 
oo 42695 9.86960 10.84276 


The error in the smallest root at n = 5 is less than 2/3 
of the error in the corresponding root at n= 3. 
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Let us now use the consistency determinant as an 
interpolation aid. We will use exact values for E at the 
symmetry points k = O and 1/2; though normally these values 
would be determined by using restricted basis functions. If 
we assume aone basis function form, then we have 


atb(E-k") + ick = 0 (6.46) 


where a = L-V, b = N and C = 2€ are to be determined by in- 
sisting that the solution E(k) to (6.46) pass through given 
points. Since E must be real, c = O,and it was found in general 
that dropping the Hermiticity requirement on the matrix never 
helps. If we now fit at the point E = .42695, k = O we find 


E = 42695 + k°. (6.47) 
Equation (6.47) evaluated at K = 7/2 gives E = 3.89435 instead 


of 2.41740, 
Assuming a two basis function form allows us to fit 


at the points 
E(0) = .42695 


E(7/2) = 2.46740 


3.37374. (6.48) 


E(1/2) 
The resulting dispersion relation is 


B2_(on7+,90634)E + (k'+.82339k°+.20467) = 0, 
(6.49) 
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which gives, besides the values (6.48) the root 
E(O) = .4813. (6.50) 


That is, a pair of "nested" curves similar to the result for 
n = 2 in the previous cases. 

If we now assume a three basis function form, we find 
that the dispersion relation can be made to fit at five points; 
the points used were those in (6.48), and 


ea) 
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~~ 
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9.86960 
E(O) = 10.84276. (6.51) 
The resulting dispersion relation is 
B3—(3k°+01.14315)E? + (3k '+2.32846u74115 .89971)E 
-(x°-18 .81468k "+-96.99219k7+45 .7070H) = 0, (6.52) 


The interpolation procedure works well between the given 
points and also as an extrapolation aid in the third curve, 
for example the third root of (6.52) is E(a/2) = 22.60420 as 
opposed to 22.20656; the error here is less that 1.8%. The 
Symmetry of the problem again seems to demand an odd number of 
basis functions. 

The various approximations to the three lowest curves 
of E(k) are tabulated below; - is the analytic solution, ED 


is the approximation using three polynomial basis functions, 
| is the 
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k = O and k = 7/2 on the two lowest curves, and k = O only 


on the third curve. 
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Table VI 

k Ee Ey Ey Z, 

0.0 10.84276 10.93605 10.91760 10.84276 
O.1 11.15673 11.39584  11.22067 11.17014 
O;2 11.73819 12.03927 11.79265 11.74979 
0.5.- 13.77463 14,13278  13.82293 13.80022 
a 17 .63365 18.03817 17.68934 17 68894. 
1.5 21.87569 22.44708  22.07127 22.09105 
1/2 22.20657 23.11187 22.73256 22.60420 


The Plane Hexagonal Lattice 
In this section we consider briefly the application 


of the new methods to a plane hexagonal lattice with two 
different potentials. The Wigner-Seitz cell and the Brillouin 
zone for this lattice are illustrated in Figures 1 and e 
respectively. Trigonometric basis functions for the consistency 
determinant are easy to find, since the general functions 


en-m 


sine x+my 

mr ) 

coser(2n-m x+my) Myth = Opti teint (6.53) 
V3 


have the appropriate periodicity. The second derivative of 
each of the expressions in (6.53) is (-167°/3) (m,n, ) times 
itself, where 


g(m,n) = m°—-mn-+n= (6.54) 
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We choose the basis functions to be of the form (6.53), and 
partially order them in the order of increasing g. If we 
eheore (mon).to be (0,0), (6,1), (130)5>(441)4 (132)5. (-1,2), 
(8,1) (0,2). (2,2) end (2,0) then @ 16 0,1,1,1,3,3,3,4;4 and 
4 respectively, and there are no other choices of m and n 
Which give £ < 4. The first few basis functions are chosen to 
be 


fH=e 1A? 
_— sinen(=) 


i = cos2n(=) 


3 V3 
f, = siner({ + y), 6.5 
mt sin me y) ( 5) 


and the basis functions form an orthogonal set such that 
fex?ar eas, (6.56) 


The matrix N, then, is (V3/4)I. The matrix L is diagonal, 
with 
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The matrices 


- or 

O43 - / ae (6.58) 
- of. 

Hy - ft sy (6.59) 


are zero except for elements like (2,3), (3,2), (4,5), (5,4), 
(6,7), (7,6), ete. The first few elements are 


= - 7/2 (6.60) 


and N22 = - N30 = 0 


ys = sy = - ™ 3/2 


+ 


N67 = 7 Teg = + 3/2 (6.61) 


The potential V(x,y) was chosen to be a two dimensional 
analogue of the Kronig-Penny potential. V is zero inside 
every Wigner-Seitz cell, and infinite on the boundary of the 
cells in such a fashion that the integral of V over any area 
is the length of the boundary included in the area. The 
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matrix V, then, 
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is given by 


Vi; - [t,V8,0 = 2 $ £40 4ds (6.62) 


where s is the edge of the Wigner-Seitz cell in Figure 1. 


The matrix V is clearly symmetric, end the first few elements 


are 
Vii = 3/2 
Vio = 0 
a See: (6.63) 
The consistency determinant was set up forn=H1, 2, 3, ... 7; 


and the successive values of the lowest energy for k = O are 


TABLE VII 
n E(0) 


2 00000 
2.00000 
1.97423 
1.97423 
1.94833 
1.94833 
1.92229 
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One expects the eigenfunction corresponding to the lowest 


eigenvalue for k = 0 to be even, and since fn» fy and P6 are 


sines, 


it is not surprising that they produce no improvement 
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in E,- For n = 7 the remaining six values of E vary between 
52,7 > 2s 57,2 (6.64) 


and represent very poor approximations to the second and third 
values of E(0). 

Let us now use polynomial basis functions satisfying 
periodic boundary conditions. The first three such polynomials 
may be chosen as 


eet 
f, = 20x"y> + hy? = 5x-y : by? 7 oe 
fs = (ay? 4Boaty? + Be? = Gx 


m 


~ 30x°y* - 15y' + 6x + 6y* (6.65) 


We note here that P3 is by chance invariant under all of the 
group operations which leave the hexagon invariant, and also 
that of/dn = © on the boundary of the hexagon. P35 then, is 
a "restricted" polynomial for the identity representation for 
r= "O, 

Using the three basis functions in (6.65) leads to 
a dispersion relation which is of degree 3 in E, 6 in k, and 
ky and 3 in the strength of the "delta function" potential 
If the potential. 1s set equal to zero, 


Bak*4k* (6.66) 


is one of the solutions. 
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At k = 0 the three roots for E(V = 0) are 0, 8.62, 

and 53.56. The value 8.62 is clearly not meaningful, and 

53.56. > 161°/3 = 5Ps637 ahs (6.67) 
which is the second (first nonzero) energy for k = 0, The 
reasonably good approximation to 161°/3 must be due to Ps 
being a properly restricted function. 

If the strength of the delta function potential is 
set to unity, then the three roots E(0) are E = 1.94559 (see 
Table VII), E = 15.5 (as opposed to E = 8.62 for V = 0) and 
53.59. If E is set equal to 8, the (lowest) curve of constant 
energy is very nearly given by 


Ke Ke = 8 - 1.94559 = 6.05441 (6.68) 
Typical points are 


TABLE VIII 
Kx Ky ke + ks 


Se ee een Sa 
+2 .46057 O 6.05441 


+2.13129 +1.23050 6.05652 
0 +2.46251 6.06396 


In the second row in Table VIII, + k, = N3ky, and so if the 
curve of constant energy had the symmetry of the lattice, 
then the two last entries in the colum ee ao LF would be 
equal; this is evidently not the case for this approximation. 
Let us now consider the method of restricted basis 
functions, first using trigonometric functions. For k = 0 
and the identity representation we must insist that the functions 


be invariant under all of the operations which leave the 
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hexagon invariant, and that the normal derivative of the 
functions vanish at the edge of the Wigner-Seitz cell. For 
the trigonometric functions, the boundary conditions are 
automatically satisfied once the symmetry conditions are 
applied. The functions are linear combinations of the basis 
functions given by equation (6.55), that is, of expressions 
of the form (6.53). Typical basis functions used were: 


4) = 2 
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Ps - 2(2/3)/2[ cos an(ae + ic) 4 aoe an(rs - y) 


: z 2(1/3) [cos on(Fz + 3y) + SOS 2n(as ss 3y) 
+ cos br ( 5% i y) + Gos bn 5 ck y) 


+ cos on( ae + y) + cos on( 7% - y)] (6.69) 


In general, six terms of the form (6.53) are required, but 
these are not necessarily distinct. 

The matrices H = -L + V and N were evaluated for 
n=1, 2, ... 6, and the eigenvalues in Table IX calculated. 


TABLE IX 


n=1 n=2 n=3 n=4 n=5 n=6 


E, 2.00000 1.91562 1.89980 1.89482 1.88492 1.388486 
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The higher roots are probably quite good since for large 
values of E, the problem approximates an empty lattice, and 
the basis functions used are empty lattice eigenfunctions. 

We now consider the method of restricted basis 
functions, this time using polynomials to construct the basis 
functions. Again we treat the case k = O and the identity 
representation. We are looking, then, for polynomials which 
are invariant under all of the operations of the group of 
the hexagon, and which have in addition the property that 
their normal derivatives vanish on the edge of the hexagon. 


The two polynomials 


p=r=x° +y (6.70) 


and 


6 


Q =r cos 60 42 4 6 


s = 15x 7 + 15x“y = ¥ (6.71) 


already satisfy the hexagonal symmetry conditions. The 
functions 


f= 
f, = 6p - 15 p- + 10p> + 29 

f. =p- Tp> ss Tp + 2pQ 

f, = 5p° - lop? - 15p' + 26p> + 10p*Q (6.72) 


satisfy both the symmetry and boundary conditions. We note 
here that f, in (6.72) above is identical to f of equation 
(6.65). The matrices N and H = -L + V were calculated, and 


the energies calculated for n=1, 2, 3, and 4, 
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TABLE X 
n=1 nee n=3 n=44 
EL 1.95358 1.91543 1489157 1.88682 
Ey -- 55 ..25eU5 54.74699 54.74006 
E -- -- 187 .40751 187.34800 
Ey, -- -- — 236.69607 


The lowest root ES must be greater than the true 
value since a variational procedure was used, and the poly- 
nomial result 1.88682 at n = 4 is better than the trigonometric 
result 1.89482, also at n = 4. The roots E, and E. from 
Table IX should be much better than the corresponding roots 
in Table X, and the apparent convergence of these roots in 
Table X must be due to very small components in the third and 
fourth polynomials of the second and third eigenfunctions of 
the problem. 

Let us now consider the plane hexagon with a potential 
which is perhaps a better two dimensional analogue of a real 


erystal. The potential used was 


ee (6.74) 


The "4" in the exponent was chosen to make the potential at 
the edge of a cell approximately 1/10 of the potential half 
way from the center to the edge of a cell. Only nearest 
neighbours were used without, then, introducing any serious 
errors. 

The polynomials for k = O for the identity representa- 
tion given in (6.73) were used as basis functions. The 
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necessary integrals were evaluated numerically using a Gaussian 
technique and are accurate to five figures. The energies 
found are given in the following table. 


TABLE XI 
1 o 3 4 
EL -1.80584 -1.89283 =-1,91442 a1 91717 
E,, saat +48 £11555 +47 48769 +47.46620 
E. _ 25 +180.60796 +179.75436 
Ei = -- —_ +357 .67424 


The Face Centered Cubic Lattice 


In choosing a physically interesting application the 
cubic lattices, having the highest symmetry among the space 
lattices, are clearly the logical choices. Of the three 
cubic lattices, the simple cubic is discarded since no 
materials are known which have this structure. The Wigner- 
Seitz cell for the body centered cubic lattice has two types 
of faces: squares and hexagons, and for this reason is some- 
what more difficult an application than the face centered 
cubic lattice for which all faces of the Wigner-Seitz cell 
are equivalent rhombuses. 

Since monovalent metals have relatively simple band 
structures compared with divalent metals, and because of the 
attention devoted to copper in the literature, it was decided 
to apply some of the new methods to copper. (See, for example, 
Howarth (1953), Fukuchi (1956), Segall (1961), Burdick (1961).) 
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The pioneer work of Seitz (1936) and Bouckaert, 
Smoluchowski and Wigner (1936) on the symmetry properties 
of electronic wave functions in crystals allows us to find 
linear combinations of powers of x, y and z having the 
symmetries appropriate to any value of the wave vector k. 
Knowing the group of operations under which a given k is 
either invariant or transforms into an equivalent k, and 
applying the Bloch boundary conditions results usually in a 
more stringent set of boundary conditions. For copper, a 
face centered cubic crystal, the irreducible representations 
corresponding to physically interesting energy values (Segall 
(1961)) for the three most symmetric points of the Brillouin 
zone are listed in Table XII. The corresponding symmetrized 
polynomials are also listed; these correspond to the lattice 
harmonics of Von der Lage and Bethe (1947) and Bell (1954). 
Finally, for the usual Wigner-Seitz unit cell bounded by the 
planes 
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are listed the boundary conditions which a wave function y, 

of the Bloch form, must satisfy. The boundary conditions are 
those of Howarth (1953) with some additions, and the error 
for X,_ has been corrected. The notation used is that of 

Jones (1960), and of Bouckaert, Smoluchowski and Wigner (1936). 
We note that since the original problem is Hermitian, the 
diverse, more stringent boundary conditions which result from 
the application of the symmetry conditions do result in the 
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vanishing of the necessary surface integrals to allow the use 
of the variational procedure. 

A start has been made in applying the method of 
restricted basis functions to this problem. For k = 0 and 
the identity representation, no polynomials exist of degree 
less than 12 (apart trom 4 constant) whieh satiety ell or the 
conditions. By taking appropriately symmetrized linear combina- 
tions of terms like xiyJzk of degree up to and including 12, 
the application of the boundary conditions resulted in 22 
homogeneous linear equations in 21 unknowns. The equations 
were solved on the LGP-30 at the University of Alberta Comput- 
ing Center using an excellent program written by Mr. A. Cseuz, 
and modified by Mr. D. L. Hunter. The coefficients Ajj, in 
Table XIII are the solutions to these equations, and are the 
coefficients of 54 4jke where 


Besuh= xt(yJzk + gdyk) 4 yi (xdzk + 24 


*) 
ijk 
+ zt (xsy* < yox*), i os J se k 


tod 4 yiatxd -} gixtyd a 


aa = xty 
ae (6.76) 


Using symmetrized polynomials up to and including 
degree 14 resulted in 28 homogeneous equations in 30 unknowns. 
One solution to these equations will, of course, be the degree 
12 polynomial already found. As yet the equations have not 
been solved, since the capacity of the LGP-30 has been exceeded. 
It was found that conventional floating point arithmetic methods 
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TABLE XIII 


~1,756,016,640 
2,004, 355,584 
399,877 632 
-43, 464, 960 
43815870, 720 
2, 347, LO7 ,O40 
3,204,432 
24,992,352 
68,690,160 
409,812,480 
-155,, 232 

shoe 20 
-3,880,800 
-23, 284,800 
-57 ,O47, 760 
3,520 

OQ 

Tes eo 
436,590 
126,126 
1,891,890 

2 401,245 
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resulted in far too little accuracy, and the program used does 
only integer arithmetic and gives exact answers. The 28 x 30 
set of equations results in large enough numbers to cause 
overflow when single precision arithmetic is used, and the 
memory of the machine is exceeded when one attempts to use 
double precision arithmetic. 

As can be seen from Table XII the determination of 
restricted polynomials, except for those corresponding to 
k = O and the identity representation, requires the satis- 
faction of boundary conditions on at least two nonequivalent 
faces. At least one of these boundary conditions always 
turns out to involve the simultaneous symmetry of W and the 
antisymmetry of ow /on, or vice versa, about a line of symmetry 
on the face of the Wigner-Seitz cell. This situation results 
in a very large number of simultaneous equations. Sets of 
such equations for several of the irreducible representations 
listed have been found. As in the case of the L polynomials, 
overflow was eventually encountered before complete solutions 
were obtained. It was found, however, that there was a 
considerable amount of linear dependence among the sets of 
equations. 

We hope to solve all of these sets of equations when 
the magnetic tape units currently being installed are fully 
integrated into the IBM 1620 system at the Computing Center. 

Another approach to the problem of constructing 
polynomial basis functions is to choose for the application 
of the boundary conditions, instead of the usual Wigner-Seitz 
rhombohedral dodecahedron, a cube of four times the volume 
bounded by the faces x = + a/2, y=+a/2, z=+a/2. Then 
at least for all of the irreducible representations listed 
in Table XII, the boundary conditions are the simple ones W = 0 
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or oy/dn = O. In compensation the cubic cell would contain 
not only the ion at the center of the cell, but also an ion 
at the midpoint of each edge. 

By using the restricted basis functions consisting 
of both polynomials and trigonometric functions, we should be 
able to determine whether or not the polynomials are more 
suitable than the trigonometric basis, as was originally 
thought probable when the investigation began. And by using 
the resulting values of E(k) for the centers of all faces 
in the first Brillouin zone, the determinant interpolation 
procedure will provide a properly symmetrized dispersion 
relation. 

Another possible interpolation procedure consists 
of fitting the special values of E(k) to a linear combination 
of polynomials in k, each having complete cubic symmetry and 
vanishing normal derivative on the surface of the Brillouin 
zone. 

Some of the above approaches may warrent further 
investigation. 
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CHAPTER VII 


SUMMARY AND CONCLUSIONS 


The new methods developed in Chapter V seem to have 
some advantages over the standard methods presented in 
Chapter IV. The standard cellular methods all make spherical 
approximations, either to the potential inside part of the 
Wigner-Seitz cell, or to the unit cell itself, or both. The 
new methods obviate such spherical approximations which, in 
some cases at least, may result in relatively large inaccura- 
cies. The Tight Binding and Orthogonalized Atomic Orbitals 
methods are essentially variational or equivalent matrix 
approaches with particular choices for the basis functions: 
atomic orbitals or linear combinations of atomic orbitals 
chosen to make the overlap integrals vanish. These choices 
of basis functions may not always be realistic. 

The plane wave method, also a variational method, 
has the obvious difficulty of apparent convergence to a 
valence state, followed eventually by convergence to the inner- 
most core state, if the potential used does not take into 
account the tightly bound core electrons. The orthogonalized 
plane wave method developed by Herring (1940) does not really 
overcome the difficulty, it merely postpones it. If the core 
states to which the plane waves are made orthogonal were 
exactly eigenfunctions of the Hamiltonian, then the variational 
procedure would indeed converge to the lowest valence state. 
Since the core states available are only approximations to the 
lower eigenfunctions of the Hamiltonian, then orthogonalizing 
the basis functions to these core states provides one with a 
set of basis functions all with relatively small, but nonzero, 
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components of the lower eigenfunctions of the Hamiltonian 
being used. The variational procedure, then, will still 
eventually provide the innermost core electron energy, but 
many basis functions will be required before this occurs. 

The convergence in the energy to a valence state with only a 
few basis functions is still only an apparent convergence, a 
"plateau" in the decrease of the lowest eigenvalue with 
increasing n, the number of basis functions. The difference 
between this plateau and the first valence state of the 
Hamiltonian depends, of course, on the core states used in 
the orthogonalization procedure and the exact core states of 
the Hamiltonian under consideration. As the core states being 
used approach the true core states, the difference between 
the plateau in the eigenvalue versus n curve and the first 
valence state will decrease, and at the same time the plateau 
Will lengthen and become an asymptote as n-o. 

For a given set of approximate core wave functions, 
there is, of course, a "best" n, for which the lowest eigen- 
value will most closely approximate the first valence state, 
but attempting to determine this best n would not appear to 
be a fruitful approach. Since the states favored by the 
variational procedure depend only on the basis functions 
used, in principle one could find the first valence state, 
say the m'th state, by using basis functions containing 
sianifiéant components of ail of the first n States, n > mM. 
One would then use the variational procedure, and pick the 
m'th lowest eigenvalue; the first n atomic orbitals might 
serve nicely as such basis functions. The matrices involved, 
of course, would be fairly large, with corresponding computa- 
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Another method of achieving the first valence state 
is by the use of a pseudopotential (Calloway (1958)) which 
does take into account the core electrons. This procedure 
would involve considerably less labour, but the accuracy of 
the resulting eigenvalues would, of course, depend on the 
accuracy of the pseudopotential. 

The three new methods developed in Chapter V appear 
to have some advantages over the standard procedures. The 
restricted basis functions, once found, may be used for the 
same lattice structure with a variety of potentials. The 
convergence of the lowest state is fairly rapid in the cases 
considered in Chapter V1, whether the basis functions are 
trigonometric or polynomial, and one cannot on this evidence 
choose between them. The trigonometric basis functions are 
certainly easier to find, though it may be that for a real 
crystal one set of basis functions may be significantly better 
than the other for rapidity of convergence. 

The consistency determinant provides a straightforward 
procedure for obtaining an approximation to the entire disper- 
sion relation. The relation obtained, however, is implicit, 
and it is possible that an inordinate number of basis functions 


would have to be used to obtain sufficient accuracy. Better 
results for the lowest eigenvalue for the plane hexagon at 

k = O were obtained with two restricted polynomials (1.91543) 
and with two restricted trigonometric function (1.91562), than 
with seven basis functions for the consistency determinant 


(1.92219). (The best answer obtained was with six trigonome 


: a O00} a dy i 3 7 4 . 1 : 1 s 4 
functions: 1.88486) It seems unlikely that the correc 


energy is less than 1.003, in which case the result fron 
consistency determinant is in error by less than 2.1%, a 
could find the entire lowest sheet of E = E(k), 


much the same accuracy. 
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From the evidence in Tables IV, V and VI, the consis- 
tency determinant used as an interpolation aid appears very 
encouraging. Together with the method of restricted basis 
functions for finding values of E(k) at symmetry points in the 


first Brillouin zone, the consistency determinant may prove 


useful. 
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